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quadrilateral.  Therefore,  the  first  part  of  this  research  program 
has  focused  on  developing  an  understanding  of  the  mechanisms  of 
rank  deficiency  and  locking  and  of  providing  a  rational  and  effective 
method  for  the  control  of  kinematic  modes. 

_ The  major  outcome  of  this  research  has  been  the  development  of 
an  effective  procedure  for  controlling  the  kinematic  modes  in  the  four 
node  quadrilateral  plate-shell  element  when  one  quadrature  point  is 
used.  In  addition,  the  insights  gained  from  that  work  have  enabled 
us  to  develop  a  triangular  plate  shell  element  with  three  nodes  which 
only  requires  a  single  quadrature  point.  Both  of  these  results  make 
possible  highly  efficient  nonlinear  transient  calculations  for  they 
permit  the  use  of  very  simple  elements  without  any  deleterious, 
effects  on  the  rate  of  convergence. 

In  addition,  some  unusual  behavior  of  higher  order  el^nents, 
such  as  the  9  node  plate-shell  element,  which  we  call  membrane  locking, 
has  been  discovered  and  investigated.  This  is  a  phenomenon  asso¬ 
ciated  with  curved  elements  that  can  lead  to  severe  errors  if  the  number 
of  quadrature  points  is  too  high.  The  importance  of  this  finding  is  that 
many  finite  element  workers  recommend  using  more  quadrature  points  when 
the  material  is  nonlinear,  in  order  to  represent  the  material  nonlinear¬ 
ity  effectively.  These  findings  show  that  if  recourse  is  taken  to 
higher  order  quadrature  in  such  elements,  the  performance  of  the 
element  may  in  fact  deteriorate  because  of  the  onset  of  membrane 
locking. 

The  results  of  the  research  conducted  so  far  indicate  that  in  the 
analysis  of  curved  shells,  an  optimal  integration  scheme  is  associated 
with  each  element  and  deviations  from  this  optimal  integration 
scheme  can  lead  to  significant  errors.  This  is  of  considerable  impor¬ 
tance  in  the  use  of  curved  shell  elements  in  structural  analysis 
where  closed  form  solutions  are  often  not  available,  because  the 
damaging  effects  of  over  and  underintegration  often  are  not  readily 
apparent.  Results  of  the  research  conducted  so  far  indicate  that 
optimal  integration  schemes  in  the  four  node  and  nine  node  elements 
are  all  associated  with  kinematic  modes.  An  effective  method  for  the 
control  of  these  modes  for  linear  problems  has  been  developed  for  the 
four-node  element.  The  performance  of  the  hourglass  control  method  has 
been  examined  in  linear  and  nonlinear  problems.  Although  some  unresolved 
difficulties  remain  to  be  dealt  with  in  nonlinear  material  problems, 
the  procedure  is  quite  effective  and  yields  a  highly  efficient  element 
which  is  suitable  for  many  applications  in  transient  analysis. 
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ABSTRACT 


A  finite  element  formulation  and  algorithm  for  the  nonlinear  analysis  of 
the  large  deflection,  materially  nonlinear  response  of  impulsively  loaded 
shells  Is  presented.  A  unique  feature  of  this  algorithm  is  the  use  of  a 
bilinear  four  node  quadrilateral  element  with  single  point  quadrature  and  a 

O 

simple  hourglass  control  which  Is  orthogonal  to  rigid  body  modes  on  an  element 
level  and  does  not  compromise  the  consistency  of  the  equations.  The  geometric 
nonlinearities  are  treated  by  using  a  corotational  description  wherein  a 
coordinate  system  that  rotates  with  the  material  is  embedded  at  the 
integration  point;  thus  the  algorithm  is  directly  applicable  to  anisotropic 
materials  without  any  corrections  for  frame  invariance  of  material  property 
tensors.  This  algorithm  can  treat  about  200  element -time-steps  per  CPU  second 
on  a  CYBER  170/730  computer  in  the  explicit  time  integration  mode.  Numerous 
results  are  presented  for  both  elastic  and  elastic-plastic  problems  with  large 

4 

strains  that  show  that  the  method  in  most  cases  is  comparable  in  accuracy  with 
an  earlier  version  of  this  algorithm  employing  a  cubic  triangular  plate-shell 
element,  but  substantially  faster. 


I .  INTRODUCTION 


In  the  nonlinear  analysis  of  impulsively  loaded  shells,  algorithms 
employing  explicit  time  integration  offer  significant  advantages  in  both  the 
economy  and  reliability  of  computations.  However,  it  is  very  important  in 
these  algorithms  that  the  bending  elment  be  quite  efficient,  because  thousands 
of  element  nodal  force  evaluations  are  typically  required  in  a  computation. 

In  an  earlier  work,  Belytschko  and  Marchertas  [1]  developed  an  explicit 
shell  code  SADCAT  based  on  the  Bazeley,  et  al  element  [2].  Although  this 
element  is  nonconforming  and  meets  the  patch  test  only  for  restricted  element 
arrangements  [3],  the  code  proved  quite  successful  and  economical  in  many 
applications.  In  addition  to  explicit  time  integration,  part  of  the 
efficiency  could  be  ascribed  to  a  corotational  formulation  which  considerably 
simplified  the  basic  equations  on  an  element  level  as  compared  to  the 
Lagranglan  formulations  which  were  then  popular.  In  Ref.  [1],  a  computational 
algorithm  was  first  developed  which  employed  vectors  to  track  the  rotations  of 
nodes  and  elements  so  that  the  arbitrarily  large  rotations  could  be  treated; 
this  was  also  applied  to  frames  In  [4], 

However,  an  element  which  does  not  meet  the  patch  test  for  all 
configurations  is  inherently  unacceptable  in  a  general  analysis  program  so 
numerous  other  elements  have  been  tried.  Experiments  with  the  Razzaque-Irons 
[5]  element  showed  it  was  too  expensive  for  explicit  computations.  Attempts 
with  higher  order  quadrilateral  elements  proved  equally  disappointing. 

In  this  paper  we  will  report  on  the  application  of  the  bilinear,  4  node 
quadrilateral  shell  element  with  one  point  quadrature,  as  proposed  by  Hughes, 


et  al  [6]  under  the  name  Ul,  which  appears  to  have  the  necessary  ingredients 
of  simplicity,  versatility,  and  reasonable  accuracy.  In  the  context  of  an 
explicit  time  integration  code,  a  simple  element  appears  to  best  because  it 
provides  the  most  accuracy  for  a  given  amount  of  computer  time;  higher  order 
elements,  while  more  accurate  for  a  given  mesh,  contain  very  high  element 
frequencies  which  severely  limit  the  stable  time  step.  Furthermore,  they  add 
substantially  to  the  complexity  of  the  computational  scheme.  This  statement 
will  be  partially  substantiated  by  comparing  the  Ul  element  to  triangular 
elements  with  cubic  fields;  other  studies  are  underway. 

Even  with  simple  elements,  reduced  integration  is  imperative  in  an 
explicit  time  integration  code  because  a  large  part  of  the  computational  cost 
arises  from  evaluating  the  constitutive  law  at  the  integration  points.  Thus 
we  have  found  that  going  from  a  single  point  quadrature  to  a  2  x  2 
reduced/selective  quadrature  more  than  doubles  the  running  time  of  the  program 
for  an  elastic-plastic  material. 

Unfortunately,  both  reduced  integration  and  selective/reduced  integration 
in  a  bilinear  plate  element  permit  zero-energy  or  kinematic  modes,  as  seen 
from  [6]  and  [7],  These  modes,  which  are  called  hourglassing  in  the  finite 
difference  literature,  will  often  quickly  destroy  a  solution.  We  will  here 
describe  the  application  of  an  hourglass  control  based  on  the  work  of  Flanagan 
and  Belytschko  [8].  The  essential  feature  of  this  hourglass  control  Is  that 
It  is  orthogonal  to  the  straining  and  rigid  body  modes  on  an  element  level 
similar  to  the  stabilization  matrix  scheme  proposed  by  Belytschko,  et  al  [9] 
for  the  selective  reduced  integration  element.  Hence  its  effects  on  the 
solution  is  minimal,  although,  in  spite  of  Its  orthogonal ity  on  the  element 
level,  it  does  slightly  stiffen  the  overall  response.  The  Ul  element  was 
first  used  with  hourglass  control  by  Taylor  [10],  who  employed  the  method  of 


Kosloff  and  Frazier  [11]. 

Another  feature  of  this  work  is  the  use  of  a  corotational  velocity-strain 
formulation.  A  corotational  finite  element  formulation  is  here  defined  as  any 
formulation  where  effects  of  rigid  body  rotation  of  the  material  are  treated 
by  embedding  a  coordinate  system  in  the  element  or  at  each  sampling  point  of 
the  element.  This  provides  a  simple  expedient  for  avoiding  the  complexities 
of  nonlinear  mechanics,  for  in  the  corotational  system  the  rate  forms  of  the 
kinematic  and  kinetic  relations  are  basically  linear  and  frame-invariant. 

The  attractive  simplicity  of  a  corotational  formulation  and  its  natural 
compatibility  with  the  finite  element  method  were  first  recognized  by  Argyris, 
et  al  [12],  who  cast  their  formulation  in  terms  of  the  natural  deformation 
modes  of  the  element.  Wempner  [13]  subsequently  developed  a  shell  theory  on  a 
similar  premise.  In  [14],  [15],  [16],  and  [12],  corotational  methods  were 
developed  and  applied  to  the  nonlinear  analysis  of  beams  and  shells  for  both 
static  and  transient  nonlinear  problems. 

The  formulation  in  [1]  and  [15]  is  a  corotational  stretch  formulation, 
for  the  strain  tensor  defined  there  corresponds  exactly  with  the  "right 
stretch"  tensor  commonly  used  in  nonlinear  continuum  mechanics;  see  [17]  to 
[19].  A  disadvantage  of  this  formulation  is  that  the  conjugate  stress  is  not 
the  physical  stress  (Cauchy  stress),  but  the  first  Ki rchhoff-Piola  stress. 

This  is  awkward  for  computer  software  because  this  stress  tensor  is  not 
symmetric,  and  Its  physical  Interpretation  is  not  as  clear  as  that  of  the 
Cauchy  stress.  Furthermore,  constitutive  models  today  are  generally  developed 
in  terms  of  physical  stress  and  its  conjugate  strain  rate,  the  rate  of 
deformation  or  velocity-strain,  so  it's  most  efficient  to  perform  element 
operations  in  terms  of  these  tensors;  note  that  the  use  of  a  corotational 
approach  does  not  affect  the  constitutive  equation  routines  at  all. 


For  these  reasons,  our  codes  have  recently  been  cast  in  this  format,  see 
for  example  [18],  The  computational  procedure  for  this  formulation  is 
actually  even  simpler  than  the  corotational -stretch  formulations  or  the 
natural  deformation  mode  formulations  if  all  computations  are  done  directly  in 
terms  of  the  velocities  and  rates.  The  rate  formulation  does  preclude  the  use 
of  large  time  steps,  but  this  is  not  a  drawback  in  explicit  time-integration 
codes,  where  numerical  stability  usually  limits  the  time  step  to  a  magnitude 
so  that  errors  in  the  integration  of  the  rate  equations  are  negligible.  The 
major  objection  we  have  found  to  this  formulation  is  that  unless  other 
measures  of  deformation  are  computed  concurrently,  the  program  provides  no 
valid  measure  of  deformation:  the  velocity-strain  tensor  itself  is  not 
integrable  [19]. 

A  word  is  also  in  order  about  "degenerate"  shell  elements,  as  pioneered 
by  Ahmad  and  coworkers  [20],  [21],  and  recently  implemented  for  general 
nonlinear  analysis  of  shells  by  Hughes  and  Liu  [22],  Although  these  elements 
have  a  compelling  cleanliness  and  possess  the  versatility  of  being  easily 
linked  with  continuum  elements,  this  is  achieved  at  some  cost.  Because  these 
elements  use  a  full  continuum  formulations,  they  require  the  evaluation  of  a 
full  3  dimensional  constitutive  equation  at  each  integration  point  and  the 
storage  of  the  complete  set  of  state  variables  associated  with  this  3D  law. 
Shell  formulations,  on  the  other  hand  only  require  a  plane-stress  law,  which 
effectively  halves  the  state  variables  and  computations.  Thus,  in  an  area 
where  we  are  still  "compute-bound",  in  that  the  size  of  computations  is  often 
limited  by  available  computer  resources,  the  shell  elements  still  appear  more 
attractive. 

In  Section  2  of  this  paper,  we  will  define  the  kinematic  and  kinetic 
state  variables  and  relations  of  the  Mindlin  theory  for  a  corotational 


description  of  shells.  In  Section  3,  the  finite  element  equations  are  given, 
followed  by  details  of  implementation,  including  hourglass  control.  In 
Section  5,  several  studies  of  the  performance  of  this  algorithm  are  reported. 

In  Appendix  A,  a  triangular  element  is  described  which  similarly  uses 
only  one  quadrature  point  per  element.  This  element  so  far  has  only  been 
tested  in  linear  situations,  but  its  characteristics  look  quite  promising. 

The  availability  of  a  triangular  element  in  conjunction  with  a  quadrilateral 
is  quite  useful  since  the  modeling  of  many  engineering  structures  requires 
triangles . 

In  Appendices  B  and  C,  the  suitability  of  some  higher  order  elements  to 
these  problems  is  examined.  It  is  shown  that  unless  reduced  quadrature  is 
employed  in  curved  elements,  a  phenomenon  called  "membrane"  locking  is 
encountered  which  leads  to  poor  results.  On  the  other  hand  reduced  quadrature 
in  these  elements  also  leads  to  kinematic  modes,  which  need  to  be  controlled. 


II.  NOMENCLATURE  AND  GOVERNING  EQUATIONS 


The  geometry  of  the  shell  is  defined  by  its  reference  surface,  or 
midsurface,  with  coordinates  denoted  by  xm,  ym  and  zm  and  by  its  thickness 
h.  The  velocity  of  the  midsurface  vm  is  given  by 


M 


(2.1) 


where  a  superposed  dot  denotes  a  time  derivative.  The  vectors  tangent  to  the 
midsurface  are  ej  and  e2  and  a  fiber  direction  is  defined  by  i.  The  fiber 
direction  is  initially  coincident  with  e3,  where 


e3  "  Si  x  e2 


(2.2) 


and  the  angle  between  £  and  e3  Is  assumed  to  remain  small,  so  that 

|  S3  i  '  1  |  <  *  (2-3) 

where  the  order  of  6  depends  on  the  magnitude  of  the  strains  and  the  accuracy 
expected;  for  most  elastic-plastic  engineering  calculations,  values  of  <s  on 
the  order  of  10"^  are  acceptable. 

The  triad  e^,  e2,  and  e3  will  be  defined  to  be  corotational  in  the  sense 
that  it  rotates  with  the  material  except  that  the  vectors  ^  and  e2  remain 
tangent  to  the  midsurface;  if  the  condition  (2.3)  is  met,  the  difference 


between  the  rotation  of  the  material  and  the  triad  e,  should  be  small.  The 

^  1 

location  of  e^  and  e3  in  the  midplane  will  depend  on  the  material  rotation  as 
defined  subsequently.  Whenever  the  components  of  a  tensor  are  expressed  in 
terms  of  the  base  vector  e^ ,  it  will  bear  a  superposed  "hat",  as  for  example 

A 

the  stress  a.  The  base  vectors  of  the  flobal  system  will  be  denoted 
by  ej.  e9  and  e§. 

In  the  Mindlin  [23]  theory  of  plates  and  shells,  the  velocity  of  a  point 
in  the  shell  is  defined  by  the  velocity  of  the  midsurface  vm  and  the  angular 
velocity  vector  9  by 

v  =  vm  -  ze3  x  e  (2.4) 


The  corotational  components  of  the  velocity  strain  (rate-of-deformation)  d  are 
given  by 


d 


ij 


3v, 


3vi 

+  - i.  ) 

3Xi  ' 


(2.5) 


Substituting  (2.4)  into  (2.5)  gives  the  following  equations  for  the  velocity 


The  velocity  strains  are  arranged  in  a  column  matrix 


subdivided  as  follows 


10b) 

10c) 


where  o'  are  the  inplane  stresses  and  a"  are  the  transverse  shear  stresses. 

A  A 

The  velocity-strain  dz  is  computed  from  the  assumption  that  az  =  0;  the 

A  A 

stresses  and  o  are  treated  primarily  as  penalty  parameters  to 
approximate  the  condition  (2.3)  and  are  not  necessarily °computed  by  the 
stress-strain  law  which  governs  the  in-plane  stresses.  This  simplifies  the 
structure  of  the  material  law  subroutine  with  apparently  no  loss  in 
accuracy . 

Note  that  the  stresses  are  always  computed  in  terms  of  corotational 
components  defined  by  the  base  vectors  e^ .  This  triad  rotates  exactly  with 
the  material  except  for  the  out-of-plane  rotation  due  to  the  difference 
between  the  rotation  of  e3  and  i,  which  is  assumed  to  be  small.  Thus  any 
anisotropic  material  law  can  be  expressed  in  rate  form  directly  as 

A  A  A  A 

a  -  S  {a,  d)  (2.11) 

without  any  corrections  for  frame  invariance.  By  contrast,  a  Jaumann  rate 
formulation  would  require  equivalent  correction  terms  for  S,  as  exemplified  in 
kinematic  hardening  models  used  by  Key  [24],  where  the  kinematic  hardening 
term  nust  be  updated  by  a  Jaumann  rate  identical  to  that  used  to  update  the 


stresses. 


III.  FINITE  ELEMENT  EQUATIONS 


The  finite  element  equation  of  motion  are  [15] 


M  v  =  fext  -  fint 


(3.1) 


where  M  is  the  mass  matrix,  fext  and  fint  the  nodal  force  matrices  arising 
from  external  forces  and  the  internal  element  resistances  respectively.  The 
internal  forces  are  obtained  by  a  topologically  appropriate  summation 


(3.2) 


A  0 

where  the  nodal  forces  f  and  moments  m  of  element  e  are  given  by  the 
principle  of  virtual  power 


..eT  _e  .  -..eT  fe  _  *eT  ^e  ,  .“eT  ie 
66j  mT  +  6Vj  fj  <S8j  rrij  +  &v  fj 


(3.3) 


A  T  A 

*  /  a  $d  o  dV 

V  '  ' 


where  Ve  is  the  volume  of  element  e;  repeated  upper  case  subscripts  are  summed 
over  the  nodes  of  the  element,  and  m®  and  f®  are  given  by 


(3.4a) 


(3.4b) 
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As  stated  in  the  introduction,  we  will  confine  ourselves  to  the  bilinear  4 
node  quadrilateral  element  with  single  point  quadrature;  by  single  point 
quadrature  we  refer  to  the  fact  that  only  one  evaluation  of  the  moments  and 
membrane  forces  is  made  within  an  element;  for  elastic-plastic  materials, 
several  integration  points  may  be  necessary  at  this  point  in  the  z  direction 
to  evaluate  the  moment  and  membrane  forces  from  the  stresses. 

The  reference  surface  of  the  shell  is  approximated  in  both  the 
underformed  and  deformed  states  by  the  elementwise  interpolation 


Nj  U.n)  t  y1 


(3.5) 


where  xj,  yj,  and  zj  are  the  coordinates  of  node  I. 
Note  that 


(3.6) 


and  the  superscript  "m"  will  be  omitted  in  the  remainder  of  this  paper  because 
all  nodal  variables  pertain  to  the  midplane 


*1  a  ^  U_C)  U~n) 


(3.7a) 


<2  (1+5)  (1-n) 


(3.7b) 


TO 


-  A 


<3  (1+5)  U+n) 


(3.7c) 


\  (1-5)  d+n) 


(3.7d) 


The  velocity  of  the  nridsurface  and  the  angular  velocity  is  approximated 
by  the  same  shape  functions,  so 


vm  ■  NjU.n)  Vj 


(3.8a) 


0  -  Nj(5,n)  8j 


(3.8b) 


Upper  case  indices  pertain  to  the  nodes  of  the  element  and  when  repeated,  as 
in  the  above,  are  summed  over  the  nodes  of  the  element.  The  velocity  strains 
at  5  *  0,  n  =  0  can  be  shown  through  Eqs.  (2.6)  and  (3.8)  to  be  given  by 


dx  *  BUvxI  *  zBlI  eyl 


“y  ■  82I  vyl  -  zB2I9xI 


2dxy  '  B2I%I  *  Bll'yl  *  z<B2I9yI  ’  Bll9xl> 


(3.9) 


Zdx2  =  Bllvzl  *  NI9yI 


2d  _  =  v_.  -  N.  e  . 


(3.10) 
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If  the  velocity  strain  matrix  is  partitioned  identically  to  the  stress  matrix, 
see  Eq.  (2.10),  then  Eq.  (3.3)  can  be  written  as 


saf  mf  +  avf  fj 


«IT«I  *  IIT  «ll 

=  /  (6d  o  +  <Sd  o  )  dV 

J  ve  '  " 


(3.11) 


where  7  is  the  shear  factor;  it  will  be  treated  as  an  arbitrary  parameter  in 
the  present  context  since  it  serves  primarily  as  a  penalization  to  enforce  the 
Kirchhoff  normality  condition  as  the  shell  becomes  thin. 

By  using  the  arbitrariness  of  the  variation  and  Eqs.  (3.9)  and  (3.11)  and 
one  point  quadrature,  we  obtain  the  following  formulas  for  the  nodal  forces 


fxj  x  A  (Bn  •/ x  +  B2I 

Kl  =  A  (B2I  /y  +  B1I  Aw* 


(3.12) 


f.r  =  A 


(B1I  /  xz  +  B2I  /yz) 


mxl  =  A  CB2I  Wy  +  B1I  **y  "  ”  K  /yz] 
myi  *  A  [-  Bn**x  -  b2i  *xy)  7/yz] 


mzl  =  0 
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where 


(3.13a) 


'  /  z 


°«e dz 


(3.13b) 


and  B^j  is  evaluated  at  the  same  point  as  in  Eq.  (3.9).  Details  of  the 
formulas  are  given  in  Appendix  A. 


IV.  IMPLEMENTATION  AND  HOURGLASS  CONTROL 
A  major  goal  in  the  programming  of  this  element  was  to  exploit  the 
simplicity  of  the  element  to  obtain  relationships  involving  few  computations 
so  that  explicit  time  integration  could  be  performed  efficiently.  Since  one 
point  quadrature  is  used  in  the  element,  hourglass  control  is  necessary.  This 
does  involve  additional  computations,  but  the  techniques  of  [8]  were  adapted 
to  the  element  so  that  the  additional  cost  is  small. 

The  element  computations  are  all  performed  in  the  corotational 

AAA 

system  (x,  y,  z).  The  function  v(x,  y,  z)  is  defined  on  the  surface  in  terms 
of  reference  parameters  £  and  n,  see  Eq.  (3.8a).  Derivatives  are  obtained 
from  the  following  matrix  equation 


where  the  comma  followed  by  a  subscript  denotes  partial  differentiation  with 
respect  to  that  subscript.  Note  that  the  third  equation  in  (4.1)  simply 
indicates  that  the  derivative  of  the  function  normal  to  the  surface  must 
vanish . 

If  Eq.  (4.1)  is  written  in  terms  of  the  corotational 
coordinates  x,  y,  z,  then  3v/3z  *  0,  so  it  follows  immediately  that  the 

A  A 

differentiation  formula  is  independent  of  3z/3£  and  3z/3n.  Thus,  the  implicit 
differentiation  formula  reduces  to 


If  we  then  use  the  identity  given  in  [8],  we  obtain 


il  IK 


A  A  A  A  A 


1  *2  ‘  *4  *3  "  *1  *4  *  H  *1  *  *3 


x4  -  x2  X1  -  x3  x2  -  x4  x3  -  X1  1=2 


(4.3) 


These  formulas  are  then  used  directly  in  the  evaluation  of  Eqs.  (3.9)  and 
(3.12),  as  given  in  Appendix  A. 

For  the  purpose  of  hourglass  control,  we  follow  [8]  and  define  the 
matrix  ^  by 


*1  ’  hI  -  (hJXoJ}  Bal 


(4.4a) 


YI  =  hI  "  B1I  +  B2I^ 


(4.4b) 


where 


hr  =  [+1,  -1,  +1,  -1] 


(4.5) 


In  the  above  equations,  the  Greek  subscripts  have  a  range  of  2  and  are  summed 

A  A  A  A 

when  repeated,  x^j  *  Xj,  Xgj  *yj.  Throughout  this  paper,  repeated  upper  case 
subscripts  are  summed  over  the  nodes  of  the  element. 

The  hourglass  generalized  strain  rates  are  obtained  by 


yl  6al 


(4.6a) 
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•B 

<3  “  n  vzl 


(4.6b) 


qM 

a 


(4.6c) 


where  the  superscripts  B  and  M  denote  hourglass  modes  associated  with  bending 
and  in-plane  (membrane)  forces,  respectively.  The  corresponding  generalized 
hourglass  stress  rates  are  given  by 


*B 

^3 


where 


Eh^A 

Clsrel5r  Bal  Bal 


C2  *  rw  JLT?~  BaI  Bal 


C3  =  rM  -V-  Bal  Bal 


(4.7) 


(4.8) 


The  constants  r0,  rw,  and  rM  are  generally  given  values  between  0.01  and 

1  *  a 

0.05.  For  elastic-plastic  materials  we  let  E  =  7  C  ,  where  C . .  are  the 

i.  ao  i  J 

constants  that  relate  the  components  of  the  in-plane  stress  tensor  by 
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i 

* 


* »  a  i 

o  =  C  d 


(4.9) 


The  nodal  forces  corresponding  to  the  hourglass  generalized  stresses  are 


N 


V 

•J 


> 

> 

A 


*3 


V 


-H  nB 

m«I  *1  Qa 

f^I  '  n  Q§  (4-W) 


where  the  total  hourglass  stresses  are  obtained  from  the  rates  as  described  in 
Appendix  A. 

An  important  aspect  of  this  hourglass  control  procedure  is  that  it  does 
not  effect  the  straining  or  rigid  boc|y  modes  for  a  Hat  element;  this  is  shown 
in  Appendix  B.  Thus  if  the  velocities  correspond  to  a  rigid  body  rotation 
about  an  arbitrary  point  or  a  rigid  body  translation,  all  of  the  generalized 
hourglass  strain  rates  vanish.  When  the  element  is  warped,  a  rigid  body 

(J 

motion  does  generate  hourglass  strain  rates  q  ;  they  may  be  almost  entirely 

a 

eliminated  by  a  procedure  described  in  Appendix  B. 

In  all  of  the  computations  reported  here,  the  central  difference  method 
was  used  for  time  integration.  A  lumped  mass  matrix  was  used  for  all 
computations.  In  addition  to  an  augmented  rotation  lumped  mass,  as  proposed 
In  [6],  a  reduced  shear  factor  <  as  proposed  in  [25]  was  used  to  reduce  the 
maximum  frequency  and  hence  increase  the  stable  time  step.  This  permits  the 
increased  rotatory  lumped  mass  to  be  scaled  so  that  the  spectral  fidelity  of 
the  finite  element  mesh  is  quite  good  over  a  large  range  of  frequencies. 


4. 
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V.  NUMERICAL  RESULTS 

Several  numerical  examples  will  be  given  to  illustrate  the  performance  of 
this  quadrilateral  element.  All  the  computations  were  carried  out  on  the  CDC 
Cyber  170/730. 

Example  1:  Impulsively  Loaded  Cantilever  Beam 

A  cantilever  beam,  shown  in  Fig.  1,  is  used  to  test  the  linear  and 
nonlinear  response  of  the  quadrilateral  plate  element.  A  uniform  pressure  is 
applied  In  the  negative  z-di recti  on  as  shown.  The  results  for  the  linear 
response  obtained  with  p  *  0.01  psi  are  summarized  in  Table  1  and  are  compared 
with  other  elements  [1]  [4],  and  the  analytic  solution  given  by  Timoshenko  and 
Goodier  [26],  The  results  for  p  =  2.85  psi,  presented  in  Table  2,  involve 
large  displacements  and  highly  nonlinear  response.  Fig.  2  illustrates  the 
time  history  for  the  tip  displacement  and  compares  it  to  the  result  obtained 
with  an  8  node,  two  dimensional  isoparametric  element  by  Shantaran  et  al  [27]. 

Example  2:  Simply  Supported  Square  Plate  Subjected  to  a  Uniform  Load 

This  problem  is  described  In  Fig.  3;  all  sides  of  the  plate  are  simply 
supported.  Due  to  the  symmetry  of  the  geometry  and  loading,  only  a  quarter  of 
the  plate  was  modelled.  The  mesh  consists  of  25  nodes  and  16  elements.  Both 
elastic  and  elastic-perfectly-plastic  materials  were  considered.  The  results 
are  compared  to  those  obtained  with  a  triangular  plate  element  [1]  and  an 
analytic  solution  by  Timoshenko  [28], 


The  elastic  and  elastic-plastic  results  are  presented  in  Tables  3  and  4, 
respectively;  time  histories  of  the  deflection  of  the  center  point  are  plotted 
in  Fig.  4  for  3  and  5  integration  points  through  the  thickness.  Note  that  the 


number  of  integration  points  through  the  thickness  used  to  evaluate  Eqs.  3.13 
makes  a  large  difference  in  the  displacement;  this  is  also  clear  from  Table  4. 

Example  3:  Impulsively  Loaded  Clamped  Beam 

A  10  in  long  aluminum  beam  clamped  at  uoth  ends  is  loaded  impulsively 
over  a  center  portion,  as  shown  in  Fig.  5.  The  material  is  elastic-perfectly- 
plastic.  Experimental  results  have  been  given  for  this  problem  by  Balmer  and 
Witmer  [29],  Fig.  6  compares  the  computed  displacement  time  history  with  the 
experimental  results. 

Example  4:  Corner  Supported  Square  Plate 

The  hourglass  modes  in  static  and  free  vibration  problems  have  been 
investigated  by  3elytschko,  Tsay  and  Liu  [9],  Here  we  will  demonstrate  the 
hourglass  modes  and  their  control  for  a  transient  problem.  We  consider  a 
square  plate  subjected  to  a  uniform  load  with  point  supports  at  the  four 
corners . 

Figure  7  shows  the  deformed  shape  without  hourglass  control.  In  order  to 
see  the  deformation,  we  have  amplified  the  results  1000  times.  This  problem 
shows  that  the  in-plane  and  w-hourglass  mode  [6,8]  produces  serious  distortion 

_3 

of  the  square  plate.  After  adding  hourglass  control  rM  =  10  ,  r  *  0.03,  the 
quadrilateral  element  gives  the  expected  deformed  shapes  as  can  be  seem  from 
Fig.  8;  the  amplification  in  this  figure  is  15000. 

Example  5:  Cylindrical  Panel 

A  120°  cylindrical  panel,  loaded  impulsively,  is  shown  in  Fig.  9.  The 
problem  is  symmetric,  so  only  half  the  panel  is  modelled.  Note  that  the  ends 
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of  the  panel  are  simply-supported  and  the  boundaries  at  the  sides  are  fixed. 
An  initial  velocity  of  5650  in/sec  normal  to  the  shell's  surface  is  applied 
over  the  region  marked  Rj. 

We  used  two  meshes  to  solve  this  problem,  varying  the  number  of 
ci rcumferential  elements  from  6  to  8  with  16  elements  along  the  length  of  the 
cylinder.  The  displacement  of  the  midpoint  along  the  crown  line  of  the 
cylinder  is  compared  with  experimental  results  [30]  and  the  triangular 
elements  results  [1]  in  Fig.  10.  The  results  obtained  for  the  6  x  16 
quadrilateral  mesh  are  not  satisfactory;  this  may  be  due  to  deficiencies  of 
one-point  quadrature  for  warped  elements.  Fig.  11  gives  the  permanent 
deformation  of  the  crown  line  of  the  panel  and  Fig.  12  gives  the  deformation 
of  a  radial  cross  section  as  compared  to  the  experimental  and  the  triangular 
element  results.  Deformed  shapes  are  given  in  Fig.  13. 

Example  6:  Spherical  Cap 

The  problem  description  and  a  top  view  of  the  mesh  are  shown  in  Fig.  14; 
fourfold  symmetry  was  used.  A  uniform  load  was  applied  over  the  cap  as 
shown.  Both  elastic  and  elastic-plastic  materials  with  the  material 
properties  given  in  Fig.  14  were  considered.  The  results  for  the  center- 
deflection  time-history  are  compared  to  the  results  obtained  by  Bathe  et  al 
[31]  using  8  node,  axisymmetric  isoparametric  elements  in  Fig.  15.  Five 
integration  points  were  used  through  the  thickness  in  the  elastic-plastic 
calculations . 

CONCLUSIONS 


A  four  node  quadrilateral  applicable  to  transient  plate  and  shell 
problems  with  material  and  geometric  nonlinearities  in  an  explicit  codes  based 
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on  the  Hughes  element  [6,7]  has  been  presented.  The  element  uses  one 
quadrature  point  in  the  plane  of  the  element  and  kinematic  modes  are 
stabilized  by  an  hourglass  control.  The  performance  of  the  element  with  one- 
point  quadrature  is  generally  quite  good,  except  when  excessive  warping  is 
encountered.  The  hourglass  control  procedure  described  here  is  easily 
implemented  and  permits  one-point  quadrature  to  be  used  regardless  of  the 
boundary  conditions,  without  mesh  instabilities.  The  use  of  one-point 
quadrature,  as  compared  to  reduced-selective  integration,  enhances  the  speed 
of  the  element  substantially;  the  element  is  also  significantly  faster  than 
the  Bazeley  et  al  [2]  element  as  used  in  [1]  with  3  quadrature  points  and,  yet 
it  possesses  comparable  accuracy. 
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APPENDIX  A 


DETAILS  OF  FINITE  ELEMENT  IMPLEMENTATION 


We  describe  here  in  detail  the  procedure  for  computing  the  internal  nodal 

forces  f  for  a  given  set  of  nodal  coordinates  and  nodal  velocities. 

-e 

Throughout  this  Appendix  a  double  numerical  subscript  indicates  a 
difference:  X32  =  X3  -  X2,  vx4i  =  vx4  -  vxi* 


Orientation  of  local  base  vectors  e^ . 

The  local  £3  vector  is  assumed  to  be  the  normal  to  the 
vectors  and  r42  as  shown  in  Fig.  A.l.  The  components  of  63  are  then  given 


y31z42  “  z31y42 
S3  =  z31x42  "  x31z42 

x31y42  "  y31x42 


(A.l) 


7  2  2  ^2 

llsll  =  (si+s2+s3^ 


(A  .2) 


Two  procedures  have  been  used  for  defining  the  x  axis.  In  the  first 

A 

procedure,  x  is  embedded  in  the  element  between  nodes  1  and  2,  side  1-2; 
however,  since  this  direction  is  not  perpendicular  to  £3  normality  is 
enforced.  This  is  quite  accurate  if  the  shear  strains  are  less  than  10%. 

A 

While  defining  x  to  join  nodes  1  and  3  would  automatically  satisfy  normality, 
the  use  of  an  axis  along  the  side  is  convenient  because  the  stresses,  which 
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A  A 

are  computed  in  the  (x,  y)  system,  are  then  more  easily  interpreted.  In 
procedure  1,  e^  is  computed  by 


X21 

r2i=  y2i  (A*3) 

Z21 

-1  =  ^21  '  ^ C21  (A*4) 

tl"l1/Ns1ll  °  (A*5) 

The  matrix  e2  is  then  obtained  by 

e2  =  e3  x  ex  (A.6) 


The  components  of  v  are  transformed  to  the  local  system  by 


where  A  is  the  matrix  of  direction  cosines  between  the  global  and  local 
system.  The  current  nodal  coordinates,  Xj  must  also  be  expressed  in  terms  of 
the  local  system  by  Eq.  (A. 7)  before  proceeding  further. 

A 

In  procedure  2,  the  coordinate  x  is  not  embedded  along  the  side  1-2; 
instead  side  1-2  is  associated  with  a  coordinate  x  and  x  rotates  with  a  spin 
as  defined  by 


3v, 


3v„ 


(A  .8) 


3x 


3y 
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/  ’  - -  -  Vi^  •  v. >. ■  y.- vfi.- . 


-> .> 


which,  by  Eqs.  (3.8),  (3.11)  and  (4.3),  gives 


A  A 


3  7K  ^24  vyl3  +  y31  vy24  '  x42  vxl3  "  x13  vx24^ 


(A.9) 


The  rate  of  the  angle  $  ,  see  Fig.  A.l,  is  then  given  by 


|  =  o)2  -  “>z(£2l)  (A. 10) 

A 

where  w  (rgj)  “*s  the  angular  velocity  of  side  1-2,  which  is 

“Z(C21 )  =  ^Vy2  ”  vyl^l  lr2lH  (A»H) 

___  ^  A  A 

The  direction  cosines  between  the  x,y  and  x,y  are  then  updated  by 


cos($n+1)  *  cos($n)  -  a$  si n ( <j>n )  +  j  a<(>2  cos ( <j>n ) 

(A. 12a) 

s1n(<frn+1)  *  sin($n)  +  a$  cos($n)  -  j  A$2  s i n ( ij>n ) 

(A. 12b) 

A<fr  =  At 

(A. 12c) 

Note  that  $  as  computed  by  Eq.  (A. 10)  is  at  time  step  n-^in  the  central 
difference  method  since  v  in  Eq.  (A.9)  is  at  time  step  n-*1^.  To  implement  this 


method,  cos  4  and  sin  $  must  be  stored  for  each  element  since  their  values  at 
the  previous  time  step,  time  step  n,  must  be  known  to  obtain  their  values  at 
n+1.  A  radial  return  is  used  to  normalize  their  values  i.e.. 


2 

s  *  (cos 

untl>  +  si 

(A. 13a) 

cos($n+* ) 

+  cos($n+1 

)/s 

(A. 13b) 

sinUn+1) 

♦  sin($n+* 

)/s 

(A. 13c) 

di recti  on 

cosines  are  then  modified  by 

S 

cos(*n+1) 

sin((j»n+1) 

0 

3 

il  - 

sin(«,n+1) 

cosUn+1) 

0 

;2t 

(A  .14) 

®3 

0 

0 

1 

ej 

where  the  far  right  hand  vector  is  that  which  appears  in  Eq.  (A. 7) 

Note  that  Eq.  (A. 7)  must  now  be  repeated  to  obtain  the  nodal  velocities 
in  the  correct  local  coordinate  system.  The  procedure  introduces  some  error 
because  Eq.  (A. 9)  does  not  use  the  correct  local  coordinate  system  to 

a 

compute  a>z ;  this  can  be  corrected  by  using  a  two  pass  procedure  or 

A 

storing  e1  for  time  step  n;  neither  alternative  appeared  to  be  worth  its 
additional  cost. 

The  area  A  is  computed  by 

A  *7  (x31y42  +  x24y13 ^  (A. 15) 

A  A 

Computation  of  strain  rates.  Once  6j  and  Vj  are  obtained  at  nodes  I,  I  =  1  t 
4,  from  the  global  components  by  the  transformation  (A.7),  the  strain  rates 
are  easily  obtained  through  Eqs.  (3.9)  and  (4.3).  The  following  formulas  are 


•*.  n.  •  ^  ^  v  \  v  \ 


K  *  7K  (y24%13  +  *13*x24> 


ay  =  7K  ^x42^yl3  +  *13^24 > 


ZKy  =  W  (x24  vxl3  +  x13vx24  +  y24vyl3  +  ^31vy24) 


Kx  =  ET  ^y24®yl3  +  y31®y24^ 


*y  *  7K  (x42exl3  +  x13ex24) 


2V  *  7K  ("x42®yl3  '  x13®x24  +  y24exl3  +  y31®y24) 


d  s  d  -  7  i c 
X  X  X 


d  3  d  -7k 

y  y  y 


H  s  H  -  7  tr 

xy  xy  xy 


(A. 16) 


(A. 17) 


(A. 18) 


The  strain  rates  nust  be  computed  at  a  set  of  quadrature  points  through  the 

A 

thickness,  -h/2  <  z  <  h/2,  if  a  plane  stress  law  is  used. 

The  generalized  hourglass  strain-rates  are  only  computed  once  in  an 
element  and  the  form  of  Eq.  (4)  is  used  directly. 


Stresses .  The  stresses  are  computed  by  a  plane-stress  constitutive  equation 
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a  s  C  d 


(A. 19) 


*  .  A 


The  transverse  she?resses  are  always  computed  by  an  elastic  law.  The 
generalized  hourgl;tress-f,ates  are  computed  by  Eq.  (4.6).  Note  that  all 
rates  are  at  time  n-tyj.  The  new  values  of  the  stresses  are  then  computed 

by 


*n+l 

o 


(A  .20) 


_h  A  h 

The  stresses  nuscomputed  at  all  integration  points  j  <  z  <  £  to 

obtain  and  e  generalized  horuglass  stresses  and  strain  rates  are 
<x$ 

computed  only  otr  each  element. 


Nodal  forces.  adal  force  contribution  from  an  element  consists  of  both 
the  nodal  forcing  from  the  physical  stresses,  Eqs.  (3.12),  and  those 
arising  from  tPralized  hourglass  stresses,  Eqs.  (4.10).  We  will  give 
the  nodal  forcessions  for  node  1: 

Kl  *  ”  XA2  ^  +  y1Q1 

fyl  “  g  (  +  >24  ^  +  ^ 

f2l-l  xz  +  x42  yz>  +  ^3  <A’21> 
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mxl  a  (x42  y  +  y24  xy^  '  4  *A  yz  +  Y1Q1 

myl  =  (  ’  y24  x  '  x42  xy J  +  4  *A  xz  +  Y1  Q2 

The  nodal  forces  and  moments  are  then  transformed  to  the  global  system  by  the 
inverse  of  Eq.  (A. 7)  which  gives 

fT  =  AY  f t  mT  =  AY  mT  (A. 22) 


APPENDIX  B 


PERFORMANCE  OF  HOURGLASS  CONTROL 
IN  RIGID  BODY  MOTION 


In  this  Appendix  it  is  shown  that  all  generalized  hourglass  strain-rates 
vanish  exactly  for  an  element  if  all  of  the  nodes  are  co-planar.  This  is 
crucial  for  the  performance  of  the  element  with  hourglass  control  because  if 
the  hourglass  strain  rates  do  not  vanish  for  rigid  body  motions,  nodal  forces 
are  generated  by  rigid  bociy  motions  via  the  generalized  hourglass  stresses. 

It  is  also  shown  that  this  condition  can  be  satisfied  when  the  element  is 
warped  and  the  nodes  are  not  co-planar  by  a  slight  modification  of  Eqs.  (4.6). 

In  this  Appendix  indicial  motion  is  used;  Greek  subscripts  have  a  range 
of  2,  Latin  lower  case  subscripts  have  a  range  of  3,  Latin  upper  case 
subscripts  a  range  of  4. 


Sj  «  Cl,  l,  1,  1] 

(B.l) 

following  equations  will  be  used 

h,s,  .  0 

(B.2) 

BaI*Sl  '  {.S 

(B  .3) 

8aISI  '  0 

(B.4) 
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where  Baj  is  defined  in  Eq.  (4.3)  and  hj  in  Eq .  (4.5);  6  is  the  Kronecker 
delta.  Eqs.  ( B .2)  to  (B.3)  can  easily  be  verified  by  simple  algebra. 

In  rigid  body  motion,  the  nodal  velocities  are  given  by 


vi  "  e1jk  "j  xk  +  v°  (B-5) 

where  e-j is  the  alternator  tensor,  is  the  angular  velocity  and  v?  the 
translational  velocity.  The  nodal  velocities  can  then  be  written  as 


"V1I  "  “j<e1j«XaI  *  e1j3‘hI  *  efjk  rkSI>  *  vtsI 


(8.6) 


where  r  is  the  vector  from  the  center  of  rotation  to  the  origin  of 
the  x  coordinate  system  and  the  assumption  that 


X3I  s  4hi 


(B.7) 


has  been  made  in  writing  the  second  term.  The  last  two  terms  will  be  omitted 
henceforth  since  by  Fqs.  (B.2)  and  (B.3)  their  inner  product  with  yj  will 
always  vanish.  The  nodal  angular  velocities  are  given  by 


0,- 


il  “i  SI 


(B.7) 


Using  Eqs.  (4.4),  (4.6a)  and  (B.7),  we  note  that 


■  ChI  -  <hJ  X8J>  S8l]  Vi  '  0 


(6.8) 


where  the  last  equality  follows  immediately  from  Eqs.  (B.2)  and  (B.4). 


Similarly,  from  Eqs.  (4.4),  (4.6b)  and  (B.6),  it  follows  that 


$3  ’  ChI  ‘  (hJ  *oJ>  Bal]  “j  e3jfl  *6I  "  °  (B'9) 

where  the  second  term  in  (B.6)  has  been  omitted  immediately  because  6333  = 
0.  The  last  equality  in  (B.9)  follows  directly  from  the  use  of  Eq.  (B.3). 
Using  a  similar  procedure  shows  that 


q”  -  4C  eeJ3  (B.10) 

Thus,  the  membrane  hourglass  strain  rates  do  not  vanish  when  t,  *  0,  i.e.  when 
the  nodes  are  not  co-planar.  However,  they  can  be  made  to  vanish 
approximately  by  letting 


qj  *  Yj  vxI  -  4?  ^  (B.lla) 

-  YX  vyI  +  4C  ;x  (B.llb) 

•M 

Equations  (B .11 )  do  not  completely  eliminate  qQ  in  rigid  body  motion  because 

•A 

the  z  coordinates  of  the  nodes  usually  do  not  satisfy  Eq.  (B.7). 


35 


vv 


t 


Fig.  1. 
Fig.  2. 
Fig.  3. 
Fig.  4 

Fig.  5. 

Fig.  6. 

Fig.  7. 
Fig.  8. 
Fig.  9. 

Fig.  10 

Fig.  11 

Fig  12. 

Fig.  13. 
Fig.  14. 
Fig.  15. 


FIGURE  CAPTIONS 

Cantilever  beam  (example  1):  problem  description  and  finite  element 
mesh . 

Tip  deflection  for  cantilever  beam  subjected  to  a  uniform  load  p  * 
2.85  psi . 

Simply  supported  square  plate  (example  2):  problem  description  and 
finite  element  mesh. 

Deflection  of  center-point  of  simply-supported  square  plate  for 
elastic  and  elastic-perfectly-plastic  materials  using  3  and  5 
integration  points  through  the  thickness  for  the  latter. 

Impulsively  loaded  clamped  beam  (example  3):  problem  description 
and  finite  element  mesh. 

Deflection  of  center  point  of  example  3  compared  to  experimental 
results  [29]. 

Deformed  mesh  of  corner  supported  plate  without  hourglass  control. 

Deformed  mesh  of  corner  supported  plate  with  hourglass  control. 

Problem  description  for  impulsively  loaded  cylindrical  panel, 
example  5. 

Displacement  time  histories  for  two  points  of  the  cylindrical  panel, 
example  5,  compared  to  experiment  [30]  and  earlier  computed  results 
[1]. 

Final  deformed  shape  of  the  panel,  example  5,  at  the  crown  line 
compared  to  experiment  [30]  and  earlier  results  [1], 

Final  deformed  shape  of  the  panel  at  the  cross-section  z  =  -6.28 
compared  to  experiment  [30]  and  earlier  results  [1]. 

Computer  plots  of  deformed  cylindrical  panel. 

Problem  description  for  spherical  cap,  example  6. 

Center  displacement  of  spherical  cap  for  elastic  and  elastic-plastic 
materials  compared  to  numerical  results  of  Bathe,  et  al  [31], 
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TABLES 


Table  1.  Parameters  and  results  for  linear  response  (p  =  0.01  psi)  of 
cantilever  beam,  example  1. 

Table  2.  Parameters  and  results  for  nonlinear  response  of  cantilever  beam, 
example  2. 

Table  3.  Parameters  and  results  for  elastic,  simply-supported  square  plate, 
example  2. 

Table  4.  Parameters  and  results  for  elastic-plastic,  simply  supported  square 
plate  example  2. 
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Thickness  t  =  1.  In. 

Young’s  modulus  E  =  12000.  pst 

Dsnslty  p  =  .1024X10“*  lb-s#c2/ln4 

Poisson's  ratio  v  **  0.2 

(a)  Problem  definition 


P 


time 

(b)  Pressure  load 


(c)  Element  mesh 

Cantilever  beam  (example  1):  problem  description  and  finite 


c 


Tip  deflection  for  cantilever  team  subjected  to  a  uniform  load  p 
2 .85  psi . 
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Length 
Thickness 
Young's  modulus 
Density 

Poisson's  ratio 
Yield  stress 
Plastic  modulus 
Pressure  load 


L  =  10.  In. 
t  =  0.5  In. 

E  a  107  psi 

p  **  2.588X10-4  lb-sec2/ln4  pel 
v  *  0.3 
0 r  *  30000.  P*I 
Ep  =  0.  psi 
P  »  300.  psi 


(a)  Problem  definition  and  element  mesh 


time 


(b)  Pressure  load 


Fig.  3. 


Simply  supported  square  plate  (example  2):  problem  description  and 
finite  element  mesn . 
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,%v  . 
— 


□  □  plate  element 
a  A  plate  element 


plaatic(3  points) 


plaatic(5  points) 


linear 

static 


elastic 


time  (sec.) 


1.00  1.20 


Deflection  of  center-point  of  simply-supported  square  plate  for 
elastic  and  elastic-perfectly-plastic  materials  using  3  and  5 
integration  points  through  the  thickness  for  the  latter. 
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-  □  plate  element 
O  experiment 


Deflection  of  center  point  of  example  3  compared  to  experimental 
results  [29]. 
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Young's  modulus  E  =  10.5 


em  description  for  impulsively  loaded  cylindrical  panel 


deflection 


o 


:=-6.28 


/^^CCQOTDo, 


z=— 9.42 


0.00 


s  □  plate  element(8Xl6) 
x  □  plate  element(8Xl6) 
a  A  plate  element (exie) 
©  experiment 


time  (msec.) 


Fig.  10  Displacement  time  histories  for  two  points  of  the  cylindrical  panel, 
example  5,  compared  to  experiment  [30]  and  earlier  computed  results 
[13. 


Final  deformed  shape  of  the  panel,  e 
compared  to  experiment  [30]  and  earl 


n  □  plate  element 
a  A  plate  element 
o  experiment 
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— i - 1  r 

1.20  1.80  2.4 

direction  (in.) 
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Final  deformed  shape  of  the  panel  at  the  cross-section 
compared  to  experiment  [30]  and  earlier  results  [1], 
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Clamped  spherical  cap  and  element  mesh 


Radius 

Thickness 

Angle 

Young's  modulus 
Density 

Poisson's  ratio 
Yield  stress 
Plastic  modulus 
Pressure  load 


r  *  22.27  In. 
t  =  0.41  in. 
a  =  28.87° 

E  a  10.5X10*  psi 

p  *=  2.45Xl0~i  lb-sec2/ln4  pel 

v  **  0.3 

<rr  »  24000.  psi 

Ep  =  0.21X10*  psi 

P  a  600.  psi 


Fig.  14.  Problem  description  for  spherical  cap,  example  6. 
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deflection 


Center  displacement  of  spherical  cap  for  elastic  and  elastic-plastic 
materials  compared  to  numerical  results  of  Bathe,  et  al  [31]. 


Table  1.  Parameters  and  results  for  linear  response  (p  =  0.01 
psi)  of  cantilever  beam,  examDle  1. 


Element 

Type 

No.  Of 
Nodes 

No.  Of 
Elements 

Time  Step 
it  (sec) 

No.  of 
Time  Steps 

Max.  Deflection 
(in) 

Period 

(msec) 

CPU  Time 
(sec) 

Euler 

Beam 

Element[4] 

6 

5 

1 .5  x 10”5 

400 

5.812 

20.03 

Triangular 

Plate 

Element[l] 

12 

20 

1.5x  10‘5 

400 

0.02408 

5.662 

126.03 

Quadri¬ 

lateral 

Plate 

Element 

12 

5 

1.5  x 10‘5 

400 

0.02454 

5.680 

25.80 

Analytic 

[27] 

0.025 

5.719 

Table  2.  Parameters  and  results  for  nonlinear  response  of 
cantilever  beam.,  example  2. 


Element 

Type 


No.  of 
Nodes 


No.  of 
Elements 


Time  Step  I  No.  of 
it  (sec)  [Time  Steps 


Max.  Deflection 

(in) 


Period 

(msec) 


CPU  Time 
(sec) 


Euler 

Beam 

Element[4] 


1.5x10' 


400 


6.321 


5.812 


20.03 


Triangular 

Plate 

E1ement[l] 


12 


20 


1.5  x  10 


-5 


400 


6.076 


5.537 


126.03 


Quadri¬ 

lateral 

Plate 

Element 


12 


1 .5  x  10 


-5 


400 


6.139 


5.640 


25.30 


2-0 

E1ement[28 


L 


0.2  x  10 


-5 


6.0 


5.600 


43* 


22 


5 


6600 


Table  3.  Parameters  and  results  for  elastic,  simply-supported 
square  plate,  example  2.  ... 


Element  No.  of  No.  of  Time  Step 
Type  Nodes  Elements  At  (sec) 


Triangular 
Plate  25 

Element[l] 


Quadri¬ 

lateral 

Plate 

Element 


32  4  x  l<f°  300 


16  6  x!0“6  200 


Max.  Oo fleet ion 
(In) 

i - 

Period 

(msec) 

0.1996 

*5  0 

Jt 

0.2001 

G.  995 

*0.2129 

1.070 

Element 

Type 


Table  4.  Parameters  and  results  for  elastic-plastic, 
sirnde  supported  square  Dlate.  example  2. 


NO.  Of 
Nodes 

No.  of 
Elements 

Time  Step 
At  (sec) 

No.  of 
Time  Steps 

Max.  Deflection 
(in. ) 

25 

32 

4x  10”6 

300 

0.3866 

25 

32 

• 

4  x  10"® 

300 

0.2478 

25 

16 

6x  10'6 

200 

0.2949 

25 

16 

6x  10*6 

2C0 

0.2511 

0.01104  143. 


0.01116  43.64 


3  layers  means  3  integration  point:  through  thickness 
5  layers  means  5  Integration  points  through  thickness 
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APPENDIX  A 


1.  INTRODUCTION 


In  the  analysis  of  nonlinear  problems,  particularly  transient  nonlinear 
problems,  computation  time  and  memory  are  often  crucial  factors.  Since 
element  stiffness  computations  are  repeated  many  times  it  is  advantageous  to 
have  efficient  and  simple  elements.  Consequently  much  research  is  aimed  at 
formulating  accurate  elements  with  these  characteri sties  [1-8]. 

In  the  analysis  of  thin  flexible  structures,  perhaps  the  most  promising 
approach  for  developing  simple  and  efficient  elements  is  that  based  on 
independent  approximations  of  the  rotations  and  displacements  combined  with  a 
reduced  order  of  shear  integration  [3-5,  9-16].  As  opposed  to  the  C1 
continuity  required  in  the  Kirchhoff  type  theory,  only  C^  continuity  need  be 
satisfied  In  this  approach.  Consequently  lower  order  shape  functions  can  be 
used  which  enhance  simplicity.  However,  the  use  of  the  low  order  shape 
functions  necessitates  reduced  integration  of  the  shear  contribution  to  the 
stiffness  matrix,  [3];  otherwise  the  elements  are  considerably  too  stiff. 
Fortunately  this  necessity  further  contributes  to  the  efficiency  of  the  C° 
element;  reducing  the  number  of  the  Integration  points  reduces  the  number  of 
computations,  and  along  with  the  simple  shape  functions,  results  in  extremely 
efficient  elements.  While  this  approach  also  has  its  drawbacks,  such  as 
possible  zero  energy  modes,  these  appear  only  for  certain  boundary  conditions 
and  then  can  be  effectively  eliminated  [17,18], 

In  some  cases  reduced  Integration  may  fall.  For  Instance  Batoz  et.al. 
[7]  examine  three  different  approaches  to  triangular  elements  and  in  their 
study  the  SRI  (selective/reduced  Integration)  triangular  element  was  found  to 
be  Ineffective.  On  the  other  hand  Hughes  and  Taylor  [19]  have  developed  a 
more  successful  one  point  quadrature  triangular  element  by  overlapping  two 
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nodes  of  the  quadrilateral  [20],  but  the  element's  performance  is  mediocre  for 
certain  element  arrangements. 

The  results  of  [7]  and  [19]  show  that  to  develop  a  successful  element 
it  is  not  sufficient  to  use  reduced  shear  Integration;  it  is  necessary  in  fact 
to  identify  those  mechamisms  which  are  associated  with  excessive  shear-strain 
energy  and  hence  shear  locking.  This  will  be  accomplished  by  fitting  a 
portion  of  the  deformation  to  an  "equivalent  Kirchhoff  mode";  this  portion  of 
deformation  will  be  called  the  bending  mode.  Although  accompanied  by  shear 
strains,  bending  modes  will  not  Involve  any  shear  strain  energy.  The 
remaining  portion  of  the  deformation  will  be  called  a  shear  mode.  The  proper 
definition  of  these  modes  is  crucial  for  the  development  of  a  successful 
element.  It  can  not  always  be  achieved  by  just  reducing  the  order  of 
numerical  Integration.  The  triangular  linear  plate  element  is  one  of  the  best 
examples  of  this  situation. 

In  this  paper  we  develop  a  new  triangular  element  with  linear  C°  fields 
which  Is  based  on  this  concept  of  decomposing  the  deformation  into  well 
defined  bending  and  shear  modes.  The  element  developed  here  shows  definite 
Improvement  compared  to  the  formulation  presented  In  [19].  We  also  identify 
the  source  of  the  difficulties  encountered  in  [7],  After  some  modifications, 
the  basic  ideas  presented  here  can  also  be  used  for  other  elements. 

Our  presentation  begins  with  general  considerations  regarding  the  shear 
and  bending  modes  of  deformation.  Section  3  contains  the  specifications  of 
the  problem  for  the  triangular  linear  plate  element.  Section  4  deals  with 
major  aspects  regarding  Implementation  while  In  Section  5  the  results  of 
numerical  applications  are  presented  including  a  discussion. 


2.  GENERAL  REMARKS 


The  main  issue  in  this  paper  is  to  determine  an  additive  decomposition  of 

the  displacements  and  rotations  of  a  plate  element  into  two  modes:  a 

bending  mode  -  associated  exclusively  with  bending  strain  energy  (regardless 

of  the  presence  of  shear  strains  in  this  mode)  and  a  shear  mode  -  which  is 

associated  only  with  the  shear  strain  energy.  Therefore  the  transverse 

deflection  w  and  the  rotations  e  ,  e  in  each  element  are  given  by 

x  y 


w  a  \  wi  Ni  *  \  (wi  +  wi )  ^ 


iy  ■  i 


xll 

f9xl 

b  j 

+ 

yiJ 

where  the  following  decomposition  is  Implied  for  the  nodal  variables 
9  ■  eb  +  8s  (lc) 

w  «  wb  +  ws  (Id) 


b  s 
w  •  w  +  w 


Here  N*  and  N^  are  the  shape  functions  for  the  displacements  and 
rotations,  respectively,  and  where  superscripts  b  and  s  designate  the  bending 
mode  and  the  shear  mode,  respectively.  Each  of  these  modes  can  Include  an 
arbitrary  amount  of  rigid  body  motion. 

The  decomposition  Is  chosen  so  that  the  element  behaves  as  closely  as 
possible  to  a  Klrchhoff  C*  element  in  the  thln-structure  limit.  The 
mechanical  reasoning  used  to  accomplish  this  task  will  be  presented  in  the 
next  Section.  In  the  decomposition,  the  shear  and  bending  nodal  variables  are 


V-yu . 
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linearly  related  to  the  total  nodal  variables,  so 


where  P*5,  Ps  are  linear  operators  emerging  from  the  decomposition.  They 
also  be  viewed  as  nonorthogonal  projection  operators. 

Assuming  the  sign  convention  shown  In  Fig.  1,  the  kinematlcal 
relationships  are 


Their  discretization  yields: 


e 


0 

w 


(4b) 


x  ■  8 


>  4  C5Sr,  B*] 


Z) 


with  Bb,  B*.  B^  defined  by  Eqs.  (1)  and  (3).  Equations  (2)  and  (4)  can  be 
used  to  find  the  strain  fields  in  the  bending  or  shear  mode  of  deformation. 
Thus,  in  the  bending  mode 


Bs  Pb 


(5b) 


while  in  the  shear  mode 


IV] 

►  -  Bb  PS  J 

- > 

ia> 

Liii  . 

[:] 

bJ 

i  '  1 

i!  j 

(6a) 


(6b) 


As  mentioned  earlier,  only  the  bending  mode  will  contribute  to  the 
bending  strain  energy,  and  only  the  shear  mode  will  contribute  to  the  shear- 
strain  energy.  Hence,  each  mode  Is  uniquely  associated  with  one  of  the  energy 
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terms.  The  total  strain  energy  Is 


•  y  y 

u  *  2  (  i  ('b)  2b  -b  ^  +  £  {x$)  2$  dA)  (7 

r  >  T  r  - 

«I  *  /  /  (Bbb)T  0b  Bbb  dA  +  J(Bss)T  Ds  Bss  dAK~  ’ 

2  J  \  A  ''  A  ”  /^J 


where  8bb,  Bss  are  defined  In  Eqs.  (5)  and  (6),  A  is  the  area  of  the  element. 


where  E  Is  Young's  modulus,  v  is  Poisson's  ratio,  h  the  thickness  and  ip  the 
shear  correction  factor.  Eq.  (7)  leads  immediately  to  the  conclusion  that  in 
the  present  formulation,  the  element  stiffness  matrix  is: 

K  -  /  (Bbb  )T  0b  BbbdA  +  /  (Bss  )T  0s  Bss  dA  (10) 

A  '  '  A 

The  above  outline  of  the  approach  does  not  give  any  rationale  for  the 
decomposition  Into  the  bending  and  shear  modes  nor  the  decomposition.  This 
crucial  aspect  of  the  formulation  will  be  discussed  in  the  next  Section  in 
connection  with  the  analysis  of  the  linear  triangular  element. 


3.  LINEAR  TRIANGULAR  PLATE  ELEMENT 


Consider  a  triangular  plate  element  in  the  local  coordinate  frame  shown 
in  Fig.  2.  Nodal  rotations  and  displacements  form  the  following  vectors 


ST  '  V'  «,2> 9 
wT  *  [*!•  *2,  «3] 


Linear  shape  functions  will 
1 

displacements  are 


y2’  ex3’  ®y3^  (lla) 

(lib) 

# 

be  used  within  the  element,  so  the  rotations  and 


(12a) 


(12b) 


where  Lj  are  area  triangular  coordinates. 

Using  Eqs  (3)  and  (12)  leads  to  the  following  forms  of 
matrices  Bb  and  Bs 


0 


0 


(14) 


A.J- 
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The  decomposition  of  the  total  displacement  into  its  bending  and  shear 
°  inodes  will  now  be  described.  To  this  end,  we  note  that  modifying  the  shear 
strain  energy  is  an  effective  means  of  improving  the  performance  of  the  C° 
element  in  the  thin  plate  limit,  for  this  may  eliminate  the  excessive  energy 
absorption  in  shear  which  leads  to  locking.  On  the  other  hand,  modification 
of  the  bending  energy  is  undesirable  since  it  may  introduce  additional  zero 
energy  modes  (compare  [3]  and  [18]  for  instance).  Therefore  the  bending  mode 
is  chosen  so  that  the  bending  strain  energy  in  this  element  is  unchanged  by 
the  decomposition.  Since  the  bending  energy  only  depends  on  the  nodal 
rotations  in  the  C°  element,  to  accomplish  this  we  can  Immediately  establish 
the  decomposition  of  9  as  follows 


eb  -  9, 


To  determine  the  decomposition  of  the  nodal  displacements  w,  we  first 
define  the  set  of  "equivalent  Klrchhoff  configurations",  which  are 
displacement  fields  w*(x,y)  with  curvatures  equal  to  that  of  the  element. 


(17a) 


i 

i 


(17b) 


=  K  =  0 

y  y  ,y 


2  K 

a  -W  .  K  s  9  +8 
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(17c) 


Note  that  the  above  implies  that  the  bending  strain  energies  in  the  C°  element 
and  any  equivalent  Kirchhoff  configuration  are  equal. 

The  curvatures  in  this  element  are  constant,  so  the  set  of  equivalent 

o 

Kirchhoff  configurations  is  described  by 


WK  *  *  2  [x2  “  X2X  '  x3(x3-x2Jy“*  ^  ’  y3y’  xy  "  x3y  ] 


"®xx  ‘  V  +  8 


(18a) 


where 


<  -  8“  6 
•w»  ^ 


(18b) 


with  8°  given  in  Eq.  (13);  a  ,  a  ,  and  8  are  free  parameters  which  account  for 

v ^  J 

K 

rigid  btfdy  motion.  At  node  I,  the  nodal  rotations  associated  with  w  are 
given  by 


3W*  -SxsM. 


(19a) 


2xrx2 
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The  nodal  displacements  of  an  equivalent  Kirchhoff  configuration  are 


In  order  for  the  C°  element  to  perform  well  In  the  thin-structure  limit, 
its  bending  mode  should  resemble  an  equivalent  Kirchhoff  configuration  as 
closely  as  possible.  Therefore,  we  will  define  the  optimal  bending  mode  as 
that  which  minimizes  the  following  measure  of  the  difference  between  an 
equivalent  Kirchhoff  configuration  and  the  bending  modes 

f(vv  8’  wi}  a  (^b  ‘  *K)T  (^b  *  *K)  +  (“b  -  *K)T  (*b  •  *K)  (21) 

where  eb,  eK  and  wK  are  given  by  Eqs.  (16),  (19),  and  (20),  while  wb  is  to  be 
determined.  The  nodal  displacement  vector  wb  (along  with  9b)  which  minimizes 
the  function  in  Eq.  (21)  Is  called  the  optimal  bending  mode.  It  is  clear  that 
the  above  procedure  is  In  fact  a  least  square  method. 

The  minimum  of  f  Is  Independent  of  the  rigid  body  translation,  b. 

Assuming  b  *  0  the  following  steps  are  needed  to  determine  wb.  First, 
according  to  Eqs.  (18b)  and  (20) 
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R  _  f1  0  1  0  1  0l 

5  "  Lo  i  o  i  o  lj 


(23b) 


The  values  of  a  and  a  which  minimize  the  function  (21)  can  then  easily  be 

*  J 

shown  to  be 


-  R  {lR  -  A)  9 

_  v  *^0  ^  v 


where  I,  is  the  unit  matrix  of  order  6.  It  also  follows  that  at  the  minimum 
-0 
H 

of  f,  w  *  w  .  Consequently,  according  to  Eqs.  (20)  and  (23b), 


wb  ,  WK  ,  .  I  x  (I.  -  A)  s 


(25a) 


0  0  0  0 

x2  0  x2  0 

x3  y3  x3  y3 


(25b) 


with  0  assumed  to  be  zero. 

Having  specified  eb  and  wb,  the  projection  operators  Pb  and  Ps  of  Eqs. 
(2),  can  be  defined.  In  view  of  Eqs.  (16)  and  (25a) 


-  7  *  (i6  -  i) 


Moreover,  since  Eqs.  (1)  and  (2)  imply  that 


Pb  +  P3  «  I„ 

•St  ^  VJ 
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the  operator  Ps  is 


The  above  expressions  for  operators  P^  and  P$  along  with  Eqs.  (4-6)  and  (13- 
15)  enable  one  to  find  the  stiffness  matrix  defined  in  Eq.  (10). 

In  Eq.  (21)  the  function  f  utilizes  local  x  and  y  components  of  rotations 
In  both  the  bending  mode  and  equivalent  Kirchhoff  configurations.  It  Is  worth 
noting,  however,  that  the  result  of  the  minimization  is  Independent  of  the 
local  frame.  To  explain  this,  note  that 


where 


49xl  ■  9xl  -  8xI<V  ay>  <29a> 

A0yl  "  eyl  "  8yl^ax’  ay)  (29b) 


Since  at  each  node  e. 


eyl'  ®xl’ 


9^j  and  ax,  Oy  are  components  of  appropriate 


vectors  (cf.  Eq.  (19))  so  are  48xj»  49^.  Therefore,  if  the  components  are 
taken  with  respect  to  a  different  coordinate  system,  the  function  f  in  Eq. 
(30a)  does  not  change.  In  particular  the  local  x  axis  can  coincide  with  any 
side  of  the  triangle. 


It  is  clear  from  the  above,  that  if  the  nodal  rotations  and  displacements 
in  the  total  configuration  coincide  with  those  of  an  equivalent  Kirchhoff 
configuration,  the  total  configuration  and  the  bending  mode  are  the  same. 
Consequently,  the  shear-strain  energy  vanishes  for  this  element  for  any 
curvatures.  We  can  show  that  this  is  not  true  of  the  formulation  presented 
in  [7]  where  just  reduced  integration  was  employed.  To  this  end,  assume  that 
the  nodal  values  of  the  rotations  are  given  by  Eq.  (19a)  while  the  nodal 
displacements  by  Eq.  (20).  In  this  case  the  shear  strains  at  the  centroid  of 
the  C°  element  are  (compare  Eq.  (19b)) 


where  xc,  yc  are  the  coordinates  of  the  centroid.  Since  in  the  thin-structure 
limit,  the  shear  strains  go  to  zero,  for  the  SRI  element  of  Ref.  [7] 

Yx  0*  Ty  ♦  0.  Thus  Eq.  (30)  (with  *  Ty  =  °)  imposes  two  constraints  on 

*v»  which  result  in  the  excessive  stiffness  of  the  SRI  element. 


4.  IMPLEMENTATION 

Both  the  bending  and  shear  contributions  to  the  stiffness  matrix  defined 
in  Eq.  (10)  form  9x9  matrices  which  are  referred  to  all  9  degrees  of 
freedom.  We  found  it  more  convenient  to  first  formulate  a  6  x  6  stiffness 
matrix,  referred  to  a  corotational  frame  in  wh1c.\ 


w  *  0  .  (31a) 

>  % 

9  «  T  '  (31b) 

S  ^  w 


where  T  is  the  transformation  matrix  resulting  from  condition  (31a).  All  of 
the  previous  considerations  are  obviously  valid  and  can  be  specialized  to  the 
following:  the  bending  part  of  the  stiffness  matrix  is 


/  (Bb)  Db  Bb  dA 
A 


(32) 


where  Bb  Is  defined  in  Eq.  (13).  The  shear  related  stiffness  matrix  is 

Es  *  -rerkr  [  -  i>T  5T  5  -  i>  "  <33> 

with  A  and  R  defined  In  Eqs.  (23  a,b).  The  stiffness  matrix  of  Eq.  (10)  can 
then  be  obtained  as  follows 


K  »  TT  (Kb  +  Ks)  T  (34) 

A  similar  transformation  has  to  be  performed  once  again  (with  a  different  T) 
to  obtain  the  matrix  K  in  a  global  coordinate  system. 
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The  integrands  in  both  Eq.  (32)  and  Eq.  (33)  are  constants,  so  the 


evaluation  of  the  integral  is  simply  the  product  of  the  area  and  the 
Integrand.  Computationally,  it  is  equivalent  to  one-point  quadrature. 

The  formulation  presented  in  the  previous  section,  utilizing  x  and  y 
components  of  nodal  rotation.  Is  one  of  the  frame-indifferent  formulations 
implemented  herein.  It  will  be  referred  to  as  the  LSC  (least  square, 
components)  formulation.  Another,  natural  and  also  frame-indifferent 
formulation  Implemented  here  utilizes  projections  of  the  nodal  rotations  onto 
sides  of  the  triangle.  We  will  refer  to  it  as  LSP  (least  square,  projections) 
formulation.  In  this  case 


f<V  V  s*  "l>  '  (2P  -  SPK)T  (2P-  s'")  *  V-  <35> 


where. 
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and  exj,  eyj  are  components  of  a  unit  vector  ej  parallel  to  Ith  side  of  the 
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1 


triangle,  Fig.  2.  Note  that  In  contrast  to  f  In  Eq.  (21),  the  first  term  of 
Eq.  (35)  Is  not  the  length  of  the  vector. 

Minimization  of  the  function  (35)  leads  to  the  following  result 


—  R  (I  -  A)  e 

3  '  '  ~  ' 


where 


—  t  T  “1  t 

R  «  3  (R  E'  E  R‘)  R  E1  E 

'v  %  %  %  '  %  s  -V 


Consequently,  If  R  Is  replaced  with  R,  Eqs.  (32),  (33)  and  (34)  are  all 
valid.  Moreover,  the  simple  form  of  matrices  R,  Eq.  (23b),  and  E,  Eq.  (36c), 
enables  one  to  perform  a  number  of  the  multiplications  In  Eq.  (38) 
analytically 


R  ET  E  RT  -  2 


2  2  2 
exl  +  ex2  +  ex3 


exleyl  +  ex2ey2  +  ex3ey3 


exleyl  +  ex2ey2  +  ex3ey3 

2  2  2 
eyl  +  ey2  +  ey3 


R  E*  E  - 
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2  2 
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6.  NUMERICAL  RESULTS 


In  order  to  evaluate  the  performance  of  this  element  and  to  compare  it 
with  other  elements,  several  square  and  circular  plate  problems  were  solved. 
The  parameters  for  all  examples  are  given  in  Table  1.  The  term  simply 
supported  here  means  that  only  the  transverse  displacements  are  constrained. 
The  transverse  load  is  discretized  in  a  manner  consistent  with  the  internal 
force  formulation;  for  a  uniformly  distributed  load,  one  third  of  the  total 
load  is  allocated  to  each  of  the  three  nodes  of  the  element.  Only  a  quarter 
of  the  plate  is  analyzed  in  each  case  because  of  symmetry.  In  the  results 
presented  (with  the  exception  of  the  square,  corner-supported  plate),  the 
deflection  of  the  center  of  the  plate  is  normalized  with  respect  to  the 
analytic  value  based  on  Kirchhoff  theory  [21]. 

The  results  for  the  square  simply  supported  plate  (Example  1)  are  shown 
In  Table  2,  where  A,  B,  CO  refer  to  the  various  discretization  patterns 
presented  In  Figure  3.  This  element  shows  marked  improvement  over  [19]  with 
mesh  A  and  a  slight  loss  of  accuracy  with  mesh  B,  so  this  element  is  less 
orientation-dependent.  It  should  be  noted  here  that  for  the  uniform  load,  the 
consistent  load  formulation  distributes  twice  as  much  load  to  the  central  node 
for  mesh  B  as  It  does  for  mesh  A,  which  Is  significant  for  the  coarse  mesh,  N 
■  4.-  Table  3  presents  the  results  obtained  by  distributing  the  load  to  the 
nodes  by  dividing  the  plate  Into  equal  square  areas  and  allocating  the 
resulting  load  to  each  node.  For  the  cross-diagonal  mesh,  the  results  of  this 
element  and  [19]  are  comparable. 

Results  for  the  circular  plate  are  presented  In  Table  4;  the 
corresponding  meshes  are  Illustrated  in  Figure  4.  Improvement  of  about  4%  to 
10%  is  gained,  over  [19]  for  mesh-type  A,  while  the  cross  diagonal  mesh  again 
yields  results  comparable  to  [19], 


We  have  also  compared  the  triangular  element  with  the  quadrilateral  with 
one-point  quadrature  and  a  stabilization  matrix.  Example  3  in  Table  1.  In 
both  cases  the  number  of  degrees  of  freedom  is  the  same.  The  mesh  used  is  of 
the  type  B,  N  »  64,  shown  in  Figure  3  and  the  results  for  both  uniform  load 
and  a  central  point  load  are  presented  in  Table  5  and  compared  to  [18].  The 
Kirchhoff  theory  solution  for  the  uniform  load  case  is  given  in  [21],  The 
performance  of  these  2  elements  is  quite  similar. 

The  convergence  rate  for  this  triangular  element  for  the  square  (edge 
supported)  plate  and  circular  plate  are  shown  In  Figures  5  and  6, 
respectively.  In  both  cases,  the  convergence  rate  is  somewhat  greater  than 
the  expected  value  of  2.0  [22],  but  no  rigorous  estimates  of  the  convergence 
rate  are  available. 

In  this  convergence  study,  the  following  measure  of  the  error  has  been 

used 

■e.  -  (/  e2da)1/2  (41) 

a 


where  a  Is  the  area  of  the  plate. 


e  »  w 


(42) 


A  FEM 

with  w  the  analytic  solution  based  on  Kirchhoff  theory  and  w  the  finite 

a 

element  solution.  To  slmplyfy  the  computations,  w  was  evaluated  at  the  nodal 
points  only  and  then  Interpolated  by  means  of  linear  shape  functions.  Thus 
the  following  difference  has  been  actually  used  within  each  element 


(43) 
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The  mesh  parameter  p  has  been  choosen  to  be  the  length  of  the  maximum  side  of 
the  biggest  element. 
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7.  CONCLUSIONS 


An  accurate  and  simple  formulation  for  the  C^  triangular  element  with 
linear  shape  functions  has  been  developed.  The  success  of  the  method  hinges 
on  the  identification  of  the  bending  and  shear  modes  and  the  use  of  the  least 
square  method  to  properly  separate  the  two  modes.  In  comparison  with  the  four 
node  bilinear  element  with  one  point  quadrature  [18],  the  element  shows 
comparable  accuracy.  However  in  shell  problems  this  element  may  prove  more 
effective  than  the  quadrilateral  with  one  point  quadrature  because  it  can  more 
effectively  handle  a  warped  surface. 

Since  this  method  assumes  a  constant  shear  in  the  element,  one  point 
quadrature  Is  sufficient  for  exact  integration  of  the  resulting  Integrals;  in 
fact,  no  numerical  Integration  is  needed.  Since  the  shear  distribution  in 
[19]  Is  linear,  a  point  probably  exists  within  the  element  at  which  the  shear 
strains  developed  In  [19]  are  equal  to  those  defined  in  this  paper.  Thus  if 
that  point  Is  used  for  the  reduced  shear  integration,  the  two  formulations 
would  yield  equivalent  stiffness  matrices.  The  present  formulation,  however, 
does  provide  a  rationale  for  a  selection  of  the  integration  point. 

Because  the  shear  distribution  in  this  element  is  constant,  it  has  one 
zero-energy  mode:  In-plane  rotation  of  the  upper  face  of  the  element  with 
respect  to  Its  midplane.  This  zero-energy  mode,  however,  disappears  in  any 
mesh  of  two  or  more  elements.  Thus  the  present  element  can  be  safely  used  in 
all  plate  problems. 
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TABLE  1 


Parameters  for  Example  Problems 


Example  1.  Square  Plate 

Uniform  Road,  simply  supported  edqes 


Dimensions: 
Thickness: 
Younq's  modulus 
Poisson's  ratio: 


10  in  x  10  in 

0.1  in 

10.92  x  105  psi 
0.3 


Example  2.  Circular  Plate 

Uniform  load,  simply  supported  edge 


Radius: 
Thickness: 
Younq's  modulus: 
Poisson's  ratio: 


5  in 
0.1  in 

10.92  x  105  psi 
0.3 


Example  3.  Square  Plate 

Uniform  load  and  concentrated  load,  corner  supported 


Dimentions: 
Thickness: 
Younq's  modulus: 
Poisson's  ratio: 


24  x  24  in 
0.375  in  . 
43.00  x  104  psi 
0.38 


r  % 


TABLE  2 

Center  displacement  and  moment  for  simply-supported 
square  plate  subjected  to  uniform  load. 


DISPLACEMENT 

MESH 

N 

NDOF 

•'Ref.  [19]  LSC 

LSP 

MOMENT 


B 


4 

cross  -.g 
diagonal  g4 


TABLE  3 

Center  Displacement  for  simply-supported  square  plate; 

Uniform  load;  nodal  forces  proportional  to  the  area  surrounding  the  nodes 


Cross 

diagonal 


NDOF 

DISPLACEMENT 

LSC  LSP 

_ _ _ 

16 

.730 

.729 

56 

.925 

.920 

208 

.989 

.988 

16 

.770 

.766 

56 

.962 

.927 

208 

.988 

.987 

28 

.913 

.936 

104 

1.035 

.986 

400 

.997 

.998 
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t 

T; 

•< 

TABLE  4 

ft 

V'. 

W\ 

*--t 

i 

Center  displacement  and  moment  for  simply-supported 
circular  plate  subjected  to  uniform  load. 

N 

ND0F 

3 

If 

12 

48 

42 

156 

DISPLACEMENT 

MOMENT 

Ref. [19]  LSC  LSP 

Ref. [19]  LSC  LSP 

cross 
diagonal!  48 


300  .976  .996 


TABLE  5 

Center  line  displacements  for  corner  supported, 
square  plate 


Nodal 

point 


DISPLACEMENT 


KIRCHHOFF 


.10855 

.10318 

.09743 

.0917'/ 

.08678 


.11955 

.11093 

.10214 

.09368 

.08602 


.11903 

.13884 

.11647 

.13392 

.11337 

.12711 

.10843 

.11919 

.10349 

.11057 

.09742 

.10183 

.09222 

.09339 

.08689 

.08576 

FIGURE  CAPTIONS 


Fig.  1.  Sign  convention. 

Fig.  2.  Geometry  of  the  triangular  element  in  a  local 
reference  frame. 

Fig.  3.  Discretizations  of  the  square  plate. 

Fig.  4.  Discretizations  of  the  circular  plate. 

Fig.  5.  Convergence  rate  for  the  square  plate. 

Fig.  6.  Convergence  rate  for  the  circular  plate. 
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APPENDIX  B 


1.  Introduction 

Since  the  early  years  of  finite  element  development,  the  use  of  C°  finite 
elements  for  the  analysis  of  thin  flexible  structures  has  been  very  tempting, 
as  the  C*  continuity  required  by  the  Kirchhoff  theory  is  very  troublesome  and 
requires  higher  order  shape  functions  [1-4],  However  a  rapid  development  of 
the  C®  approach  started  only  with  the  application  of  reduced  integration  [5- 
12]. 

The  effectiveness  of  this  technique  arises  from  its  elimination  of  the 
excessive  shear  contribution  to  the  stiffness  of  thin  structures  whose 
response  is  usually  dominated  by  their  bending  properties.*  However,  it  was 
not  clear  how  much  of  the  shear-rel ated  stiffness  should  be  eliminated. 
Consequently  the  elements  employing  reduced  integration  were  developed  on  a 
"trial  and  error"  basis;  certain  interpolations  and  integration  schemes  were 
usually  assumed  and  their  consequences  were  examined,  [9-12].  Elements  that 
did  not  perform  "well"  were  rejected.  Although  a  similar  approach  was  used  in 
the  development  of  Kirchhoff  C1  elements,  in  this  case,  approximations 
consistent  with  the  theory  usually  give  acceptable  elements.  Although  their 
convergence  properties  and  error  characteristics  may  vary,  they  seldom  fail, 
as  for  instance  when  the  thickness  of  the  plate  decreases.  This  is  not  true 
of  C°  elements  employing  reduced  integration  (for  'instance  serendipity  plate 
elements  [9-12])  and  perhaps  for  this  reason  the  technique  is  sometimes  viewed 
more  as  a  trick  than  a  legitimate  method.  In  mixed  methods,  [17,18],  similar 
trial -and-error  procedures  have  been  used. 

The  equivalence  between  the  reduced-integration  displacement  approach  and 
well  established  mixed  methods  [19]  contributed  significantly  to  the 

*  A  more  complex  phenomena  occurs  when  a  curved  structure  is 
analyzed  [14-16], 


legitimacy  of  the  former.  However,  doubts  remained  as  to  how  to  create 
successful  elements;  trial -and-error  approaches  are  still  in  use  and 
consequently  the  resulting  elements  are  often  considered  insufficiently 
reliable.  The  zero-energy  modes  that  often  accompany  reduced  integration  can 
also  cause  severe  difficulties. 

A  deviation  from  this  path  was  originated  by  MacNeal  [20]  who,  by 
comparing  energy  terms,  attempted  to  create  C®,  low  order  elements  of  accuracy 
equivalent  to  that  of  higher  order  elements.  A  similar  approach  was  later 
used  by  Parish  [21]  to  justify  the  use  of  reduced  integration  within  9-node 
Lagrange  plate  elements. 

Another  approach  for  improving  accuracy  and  eliminating  zero-energy  modes 
was  proposed  by  Hughes  and  Tezduyar  [22]  who  developed  a  successful  four-node 
quadrilateral  plate  bending  element.  The  idea  was  similar  to  that  presented 
in  [23]  and  consisted  In  changing  the  discrete  derivative  operator  B  to 
guarantee  good  behaviour  of  the  element  in  thin  plate  limit.  After  the  change 
of  the  matrix  B  there  is  essentially  no  need  for  reducing  the  order  of  the 
integration  to  obtain  good  res?  An  identical  approach  to  the  analysis  of 

perhaps  the  simplest  plate  bending  element,  a  triangle  with  linear 
approximations  of  displacements  and  rotations,  was  presented  by  Hughes  and 
Taylor  [24].  The  results  for  this  element  depended  very  strongly  on  the  mesh 
arrangement.  Consequently  the  reduced  integration  technique  was  applied  to 
the  shear  terms  along  with  the  modified  B  matrix  to  alleviate  the  dependence 
on  the  mesh  orientation.  However  doubts  arose  about  tne  correctness  of  the 
modified  B  matrices.  This  subject  has  been  discussed  in  [25],  where  a  new 
approach  to  the  modification  of  the  B  matrix  was  proposed. 

Another  approach  to  elements  is  to  use  reduced  integration  with 
stabilization  matrices  [26,27],  which  eliminate  the  zero-energy  modes.  These 
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approaches  have  the  same  or  better  accuracy  than  [22],  but  the  selection  of 
stabilization  parameters  for  nonlinear  problems  is  an  open  question  which  is 
probably  not  trivial. 

This  study  is  aimed  at  identifying  the  factors  which  are  most  essential 
for  the  success  of  a  C®  flexible  element.  This  is  done  by  means  of  two 
elements:  the  linear  beam  element  and  the  linear  triangular  plate  element. 

The  salient  characteristic  of  the  first  is  that  all  the  formulations  discussed 
lead  to  identical  results;  this  is  not  the  case  for  the  plate  element.  We 
believe  this  enables  one  to  clearly  see  the  major  features  of  the  problem. 

The  next  Section  contains  general  remarks  concerning  C°  flexible  plate 
elements.  The  various  formulations  for  the  beam  problem  are  given  in  Section 
3,  while  the  linear  triangular  plate  element  is  discussed  in  Section  4. 
Numerical  results  and  conclusions  are  presented  in  Sections  5  and  6, 
respectively. 
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2.  General  Remarks 


constants,  x  is  the  shear  correction  factor.  The  matrices  BK,  Bc  are  the 
discrete  forms  of  the  kinematical  relationships 


where  k  ,  k  ,  k  are  curvatures,  y  ,  y  are  shear  deformations  and  q  is 

a  Jr  Ajr  X  j 

the  vector  containing  all  elemental  degrees  of  freedom.  Positive  rotations 

8  ,  0  and  displacement  w  are  shown  in  Figure  1. 
x  y 

It  is  known  that  if  low  order  shape  functions  are  used  and  if  the 
integral  in  Eq.  (lc)  Is  Integrated  exactly,  the  plate  elements  are  too  stiff, 
or  "lock",  as  their  thickness  decreases  [6].  This  is  due  to  the  fact  that  too 
much  of  the  work  performed  by  the  applied  forces  is  converted  to  shear-strain 


energy.  To  alleviate  this  problem,  the  integral  expressing  the  shear 
contribution  to  the  stiffness  matrix  is  often  evaluated  using  quadrature 
schemes  of  lower  order  than  that  required  for  exact  integration.  By  doing 
this,  the  portion  of  the  shear  strain  energy  associated  with  higher  order 
distribution  of  the  shear  strains  Is  removed  and  the  performance  of  the 
element  is,  in  general.  Improved. 

8ut  this  is  not  always  true.  There  are  cases,  like  the  triangular  plate 
element  discussed  subsequently,  that  can  not  be  treated  in  this  way. 

Moreover,  reduced  integration  often  gives  zero-energy  modes  that  are  highly 
undesirable.  Furthermore,  one  can  never  be  sure  whether  reduced  integrator 
removes  the  correct  portion  of  the  shear  strain  energy.  It  may  therefore  be 
better  to  remove  the  troublesome  terms  by  using  mechanical  arguments.  This 
reasoning  leads  to  a  different  matrix  Bg  that  replaces  B$  of  Eq.  (lc). 
Depending  on  whether  or  not  all  the  troublesome  shear  terms  have  been  removed 
while  defining  the  B$  matrix,  the  Integral  of  Eq.  (lc)  can  be  integrated 
exactly  or  under! ntegrated  [22,24].  The  point  is  to  find  a  general  method 
that  eliminates  all  the  troublesome  terms  [25].  As  pointed  out  in  [22], 
mechanical  reasoning  may  provide  a  Bs  matrix  which  is  far  more  effective  than 
reduced  Integration. 

If  a  mixed  model  is  used  to  formulate  the  element  stiffness  matrix, 
internal  forces,  work -conjugate  with  the  strains  given  In  Eq.  (4)  and  (5),  are 
Interpolated  Independently  of  displacements 


My  r  "  ~M  &M 


(7) 


-t  3t 


where  P^,  Py  describe  the  distributions  of  M  and  T  respectively. 

For  (3^  and  gy  Independent  of  each  other  (and  this  is  often  the  case  in 
practice)  the  two  components  of  the  stiffness  matrix  are  (compare  [16]): 


(8a) 

(8b) 


Moreover,  If  for  a  given  displacement  field  the  moment  distribution  is  exactly 
the  one  that  would  occur  in  the  displacement  approach  i.e.  if 


P 


M 


3b  3b 


(9) 


then  the  shear-related  stiffness  matrix  Ks  is  still  defined  by  Eq.  (8b),  but 
equation  (8a)  becomes 


3b  *  I  3b  3b  3b 

A 


(10) 


In  this  case  only  the  distribution  of  the  shear  forces  has  to  be  defined. 

This  version  of  the  mixed  formulation  will  be  used  throughout  the  paper. 

Malkus  and  Hughes  [19]  showed  that  if  the  kinematical  description  in  the 
displacement  and  mixed  approaches  is  the  same,  then  for  each  reduced 
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Integration  scheme  applied  within  a  displacement  approach,  there  exists  an 
Internal  force  distribution  In  the  mixed  approach  such  that  the  stiffness 
matrices  obtained  by  the  two  methods  are  Identical.  In  our  case,  for  a  shear 
underintegration  scheme,  there  exists  a  shear  force  distribution  Py  such  that 
the  stiffness  matrices  defined  by  Eq.  (1)  and  Eq.  (10)  are  identical.  This 
equivalence  will  be  used  later. 

In  subsequent  chapters  we  will  be  discussing  various  equivalent 

formulations.  By  equivalent  formulations  we  will  mean  those  that  result  in 

Identical  stiffness  matrices.  In  all  of  the  problems  discussed,  rigid  body 

# 

motion  will  be  eliminated  and  only  the  corotational  stiffness  matrix  will  be 
compared.  This  diminishes  the  number  of  degrees  of  freedom,  yet  preserves 
full  generality  of  the  analysis. 


3.  Linear  Beam  Element 


Consider  a  beam  of  length  L,  width  b  and  depth  h.  We  assume  linear 
distributions  of  rotations  and  displacements  which  depend  on  four  nodal 
quantities  9^,  w^,  72,  *2*  wflere  bars  denote  quantities  defined  in  the  global 
system.  The  corotating  frame  (x,y)  Is  defined  so  that  the  transverse 
displacements  vanish,  so  the  deformation  is  totally  described  by  the  rotation 
field 


6  -  fljU-O  +  6^  ,  5  ■  £ 


(Ha) 


ST  ■  [Sj.  e2] 
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The  curvature  and  shear  strains  are  given  by 

Kx a  9*x  *r  (02  •  91) 

Yx  *  e  +  w,x  "  91  <1_5)  +  0 2^ 

The  discrete  forms  of  the  above  kinematical 
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relations  are  given  by 


It,  -r  c-1*  « 
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Is  -  [1-5,  5] 


1)  Reduced  integration  displacement  formulation 


In  this  approach  the  integrals  in  Eqs.  (lb)  and  (lc)  are  evaluated  using 
5  *  as  the  integration  point.  This  is  exact  integration  for  the  first  but 
reduced  Integration  for  the  second  integral.  The  resulting  stiffness  matrices 
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i  1 )  Mixed  formulation 

Use  of  the  Eqs  (8)  will  now  be  made  with  and  Bg  defined  in  Eqs.  (13), 
P^  given  In  Eq.  (9)  and  Py  defining  a  constant  distribution  of  the  shear 
force,  e.g. 

PT  •  Cl]  (15) 


Exact  integration  of  the  resulting  formulas  yields  the  matrix  Kb  given  in  Eq. 


[• 


(14a)  and  matrix  K$  given  in  Eq.  (14b).  A  constant  shear  distribution  in  the 
mixed  method  is  therefore  equivalent  to  the  reduced  midpoint  integration  in 
the  displacement  approach. 


The  above  two  formulations  are  well-known  [6,19],  Here  we  present 
another  formulation  based  on  a  concept  developed  as  follows.  First  note  that 
in  pure  bending,  a  Ki rchhof f-type  element  in  which  the  curvature  is  constant 
undergoes  equal  but  opposing  nodal  rotations.  If,  in  the  thin  structure 
limit,  the  present  C®  element  (which  also  gives  constant  distribution  of  the 
curvature)  is  to  behave  like  its  counterpart  In  the  Kirchhoff  theory,  the 
symmetric  part  of  deformation  that  preserves  curvature  and  is  shown  in  Fig. 
2b,  should  be  associated  with  no  shear  strain  energy.  The  remaining, 
antisymmetric  deformation.  Fig.  2c,  does  not  change  the  curvature  and  is 
characterl zed  by 

i\  m  "7  (01  +  02)  <16) 

Associating  only  this  mode  with  the  shear  strain  energy  we  arrive  at  the 
matrix 


which  reflects  the  shear  strain  distribution  In  the  antisymmetric  part  of 
deformation.  Moreover,  since  the  relationship  between  Bs  of  Eq.  (17)  and 
B$  of  Eq.  (13b)  is 


(18) 


Is  ■  Is  <«  ■  i> 

use  of  B$  In  Eq.  (lc)  instead  of  Bg  and  exact  integration  is  equivalent  to 
the  reduced  midpoint  integration  presented  as  formulation  (i).  Consequently 
once  again  the  matrix  K$  is  given  in  Eq.  (14b).  The  matrix  §b  remains 
unchanged,  which  means  that  £b  is  still  given  by  Eq.  (14a). 

i v)  Displacement  formulation  with  a  modified  distribution  of  the  transverse 
displacements 

Now  the  symmetric  portion  of  deformation,  discussed  in  the  previous 
formulation  will  be  introduced  more  explicitly;  a  quadratic  transverse 
displacement  field  will  be  associated  with  it.  Thus  we  assume 

wK  -7  02  -  0j)  (1-0  5  (19) 

while  the  rotation  field  Is  still  described  by  Eq.  (lib).  This  displacement 
field  will  be  used  within  the  standard  displacement  approach,  not  employing 
any  reduced  Integration  whatsoever.  Therefore  the  problem  reduces  to  the 
evaluation  of  Bb  and  B$  and  to  the  exact  Integration  prescribed  by  Eqs.  (1). 
According  to  Eqs.  (12a, b)  <x  does  not  change  and  neither  does  Bb  while 

Yx  *  0  +  w,x  *  M1”5*  +  ®25  +7  (02  "  9l)  (1  ”  ^  (20) 

-  7  (9j  +  e2) 

This  results  in  the  following  matrix 

S-[M]  121  > 
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Since  |b  remains  unchanged  and  Bg  Is  the  same  as  §s  of  Eq.  (18),  the  last 
two  formulations  are  equivalent.  Furthermore,  as  opposed  to  the  previous 
formulation,  no  simplifications  ("tricks")  have  been  used  in  the  present  one. 


Note,  that  In  the  thin  structure  limit,  when  the  constraint  yx  *  0  is 
enforced,  Eq.  (20)  yields 

9X  -  -9 2  (22) 

This  Is  consistent  with  the  quadratic  distribution  of  the  transverse 
displacement  given  in  Eq.  (19)  and  shows  that  all  of  the  above  formulations 
attain  the  accuracy  of  the  quadratic  K1  rchhof f-type  element.  However,  the 
quadratic  Klrchhoff  element,  developed  without  any  shear  deformation,  would  be 
associated  with  a  complicated  assembly  procedure  since  the  two  nodal  rotations 
would  not  be  Independent. 

v)  Mixed  formulation  with  the  modified  distribution  of  the  transverse 
displacements 

A  quadratic  distribution  of  displacements,  Eq.  (19)  and  a  constant 
distribution  of  the  shear  deformation  can  be  used  to  show  that  the  mixed 
formulation  leads  to  precisely  the  same  results  as  those  in  the  formulation 
(1 v) .  This  can  be  Immediately  concluded  from  Fraeljs  de  Veubeke's  limitation 
principle,  [28],  since  In  the  kinematic  approach  presented  in  (iv)  the  shear 
distribution  Is  also  constant,  Eq.  (20). 

vl )  Displacement  formulation  based  on  the  optical  bending  configuration 

A  detailed  discussion  of  the  approach  is  given  in  [25];  its  basic  idea  is 


the  following.  In  order  for  a  flexural  C®  element  to  perform  well  in  the 
thin-structure  limit,  there  must  exist  a  properly  defined  deformed 
configuration  associated  with  bending  strain  energy  and  no  shear  strain 
energy.  This  configuration  will  be  called  the  bending  mode.  The  additional 
deformation  required  to  bring  this  mode  to  the  total  deformed  configuration 
will  be  called  the  shear  mode;  It  Is  associated  with  a  shear  strain  energy  and 
no  bending  energy.  Since  the  assumed  coordinate  system  is  corotational  only 
for  the  total  configuration,  the  bending  and  the  shear  mode  may  be 
characterized  by  nonzero  nodal  displacements.  However,  at  each  node,  the  sun. 
of  the  displacements  describing  the  two  modes  has  to  vanish  whereas  the  sum  c 
the  rotations  should  give  the  Initial  rotations.  Both  modes  are  described  bj 
the  linear  shape  functions.  Their  proper  definition  is  essential  for  a 
successful  development  of  a  C°  element. 

Existing  works  clearly  Indicate  that  shear  strain  energy  should  be 
modified  to  achieve  good  behavior  of  the  C°  elements  in  the  thin-structure 
limit.  Moreover,  they  Indicate  that  the  modification  of  the  bending  strain 
energy  (Introduced  for  Instance  by  reduced  bending  integration)  is  not 
desirable  since  it  usually  Introduces  additional  zero-energy  modes.  For  that 
reason,  the  bending  strain  energy  in  the  total  configuration  and  in  its 
bending  mode  should  be  the  same.  To  Insure  that  this  requirement  is 
fulfilled,  the  nodal  rotations  In  the  bending  mode  are  assumed  to  be  the  same 
as  In  the  total  configuration.  This  implies  that  the  nodal  rotations  in  the 
shear  mode  are  zero;  thus  the  shear  strains  that  should  be  taken  Into  account 
are  completely  defined  by  nodal  displacements  in  the  shear  mode. 

Even  without  a  quantitative  formulation  the  following  remarks  can  be 
made.  First,  It  Is  clear  that  the  shear  strain  described  above  involves  only 
the  first  derivatives  of  the  transverse  displacements.  The  polynomials 


describing  the  distribution  of  the  shear  strains  are  therefore  one  order  lower 
than  those  resulting  from  a  given  displacement  and  a  nonzero  rotation  field. 
Second,  although  the  first  remark  indicates  that  certain  shear  strains  have 
been  removed,  they  do  not  have  to  be  the  same  as  those  removed  by  reduced 
shear  integration.  The  present  approach  is  based  on  mechanical  reasoning  and 
can  serve  as  a  guide  to  the  appropriate  reduced  quadrature. 

To  determine  the  bending  mode,  we  consider  a  set  of  "equivalent" 

Kirchhoff  configurations:  they  are  characterized  by  a  curvature  identical 
with  that  in  the  total  configuration.  They  can  be  described  by  superposing 
the  quadratic  displacement  field  of  Eq.  (19)  and  a  linear  field  resulting  from 
a  rigid  body  motion.  Then  we  select  the  Kirchhoff  configuration  whose  nodal 
rotations  are  closest  -  in  an  average  sense  -  to  the  already  defined  rotations 
in  the  bending  mode  (being  equal  to  the  total  rotations).  By  identifying  the 
nodal  displacements  In  this  particular  equivalent  Kirchhoff  configuration  with 
the  displacements  in  the  bending  mode,  we  define  the  optimal  bending 
configuration  (or  bending  mode). 

Since  only  rotations  are  compared  in  the  evaluation  of  the  nodal 
displacements,  the  rigid  body  translation  is  irrelevant.  Thus,  rotation 
around  the  node  1  of  magnitude  a  Is  considered.  The  difference  between 
rotations  is  measured  by  a  sum  of  squares,  so  we  miminize 

♦  (a)  *  (AQ x)2  +  (A02)2  (23) 

where 
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with  respect  to  a  .  It  turns  out  for  this  case  a  can  easily  be  chosen  so 

that  4  »  0,  which  is  obviously  the  minimum 

°  *  7  (91  +  92^  (25) 

This  defines oa  particular  Kirchhoff  configuration  shown  on  Fig.  3.  To 
transform  the  configuration  to  the  original  configuration  the  nodal  points 
have  to  be  displaced  by 

AWj  *  0  Aw 2  »  aL  (26) 

These  are  nodal  displacements  in  the  shear  mode  that  result  in  the  following 
shear  deformation 

1  *  o  (27) 

In  view  of  Eq.  (25)  the  modified  matrix  B$  is 


and  It  coincides  with  the  one  given  In  Eq.  (17).  Consequently  the  present 
formulation  is  equivalent  to  all  the  previous  ones. 

We  have  Illustrated  these  concepts  In  terms  of  a  beam  element  for 
simplicity.  It  is  also  interesting  to  note  that  all  of  the  preceding  beam 


formulations  give  the  same  results.  In  the  next  Section  we  will  discuss  these 
formulations  in  application  to  the  triangular  linear  plate  element.  In  this 
case  they  are  not  all  equivalent.  The  results  will  suggest  the  formulation 
which  can  be  safely  used  in  a  wide  class  of  problems. 


4.  Triangular  linear  plate  element 

The  triangular  linear  plate  element  will  be  considered  In  the  local 
corotational  frame  of  reference.  Fig.  4,  such  that  the  nodal  transverse 
displacements  are  zero.  Its  deformation  Is  therefore  described  by  the  nodal 
rotations 
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whereas  the  rotations  within  the  element,  e  and  e  ,  are  distributed  linearly 
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with 
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and  denoting  the  area  coordinates.  Using  the  above  expressions  and 
utilizing  Eqs.  (4)  and  (5)  one  arrives  at 


Note  that  Bb  defined  in  Eq.  (32)  is  constant  while  8$  of  Eq.  (33)  is 
linear.  Thus,  one  point  centroidal  integration  gives  the  exact  value 
for  Kb  and  underintegrates  <s  .  The  expressions  defining  the  two  matrices  are 

Kk  -  A  bJ  D.  Bk  (34a) 

K  *  A  B,  ol  BCc  (34b) 


where  A  is  the  area  of  the  triangle  and  B|  is 

evaluated  at  the  centroid  (Lj  *  L2  *  t-3  ■  -). 
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the  matrix  B$,  Eq.  (33), 


1i )  Mixed  formulation 

We  assume  that  PM  is  chosen  according  to  Eq.  (9)  and  that  Kb  is  given  by 
Eq.  (10).  Both  components  of  the  shear  terms  are  independently  assumed  to  be 
constant  within  the  element,  so 


(35) 


and  then  Eq.  (8b)  reduces  to 


Since  B$  Is  a  linear  function  over  the  triangle,  one  point  centroidal 
Integration  evaluates  the  first  and  the  last  integral  exactly. 


Furthermore  D„  is  constant  so 
~s 


Ss 


A  Bf 


(A^) 


-1 


A 

-s 


A  Bf  D  B^ 
-s  ~s  ~s 


(37) 


which  is  identical  with  Eq.  (34b).  This  formulation  is  therefore  equivalent 
to  the  previous  one. 

Both  formulations  as  well  as  their  equivalence  were  reported  in  [29]. 
However  the  approach  turned  out  to  be  Ineffective.  Further  search  for  an 
appropriate  approach  to  a  linear  C°  triangular  element  has  led  to  the 
formulations  presented  In  the  following. 

Ill)  Displacement  formulations  with  the  modified  matrix  Bc 

Both  the  curvatures  and  the  shear  strains  along  side  of  the  element 
depend  on  the  projection  of  the  total  rotation  on  the  particular  side.  For 
the  linear  rotation  field  defined  In  Eq.  (31)  the  projections  0^,  i  *  1,2,3, 
are 

®1  *  ®if1(l-Ci)  +  £•)  (38) 

with 

•1.1‘iiSi  8i.hi  "sliSi  <39> 

where  for  1*3,  1+1  should  be  identified  with  node  1,  9^  is  specified  in  Eq. 


(31),  e^  are  the  unit  vectors  shown  on  Fig.  4  and  ^  (0,1)  parametlzes  the 

1-th  side  of  the  triangle.  Thus  the  curvatures  and  the  shear  strains 
associated  with  each  side  are 
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The  above  equations  are  analogous  to  Eq.  (12a, b).  One  can  therefore  draw  the 
following  conclusion:  if  the  present  plate  element  is  to  behave  well  in  the 
thin  plate  limit,  the  linear  portion  of  the  shear  deformation  for  each  side 
has  to  be  related  to  a  deformation  associated  with  no  shear  strain  energy 
(this  led  to  success  In  the  analysis  of  beams).  The  portion  of  the  shear 
strains  that  should  be  associated  with  the  strain  energy  is  therefore  (compare 
Eq.  (16)) 
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These  values,  considered  for  1  ■  1,2,3,  define  the  modified  shear  deformation 
over  the  entire  triangle 
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can  be  obtained  by  Inverting  the  relationships 

r*  -  (*?)T  Sl  -  (aUl)  Sf  ,-1-2-3  <44> 

Taking  Into  account  Eqs.{39),  (41),  (42)  and  (44)  as  well  as  the  particular 
position  of  the  reference  frame.  Fig.  4,  one  arrives  at 
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The  approach  presented  above  was  first  proposed  by  Hughes  and  Tezduyar, 
[22],  for  the  analysis  of  the  quadrilateral  element,  and  then  was  used  by 
Hughes  and  Taylor,  [24],  in  the  analysis  of  the  triangular  element.  Since 
the  shear  strain  distribution  defined  by  Ts  in  Eq.  (45)  is  linear, 
no  single-point  integration  (associated  with  constant  shear  strains  over  the 
element)  can  be  equivalent  to  the  present  approach.  Therefore  this 
formulation  is  clearly  different  from  the  two  previous  ones.  Yet,  the 
derivation  of  Bs  as  well  as  the  difference  between  Eqs.  (33)  and  (45)  clearly 
indicates  that  some  portion  of  the  shear  strain  energy  has  been  removed. 

Although  use  of  gs  instead  of  Bg  in  Eq.  (10)  should  not  require 
reduced  integration,  exact  Integration  has  been  found  to  lead  to  results  very 
strongly  dependent  on  the  mesh  orientation,  [24];  consequently  one  point 
quadrature  has  been  applied  by  Hughes  and  Taylor,  [24],  but  the  integration 
point  has  been  selected  so  that  Its  location  (and  the  results)  depend  on  the 
local  numbering  of  nodes.  Fig.  5  shows  three  locations  of  the  integration 
point  P,P',P"  for  three  different  numbering  of  nodes  123,  r2'3'  and 
1",2",3M.  It  is  therefore  clear  that  the  integration  point  has  not  been 
selected  properly.  More  Importantly,  however,  the  above  analysis  Indicates 
that  the  matrix  should  be  defined  In  a  better  way;  a  correct  definition  of 
1S  should  not  necessitate  any  reduced  integration.  These  problems  will  be 
discussed  subsequently. 

1 v)  Displacement  fromulatlon  with  a  modified  distribution  of  the 
transverse  displacements 

Here,  a  displacement  formulation  equivalent  to  the  previous  one,  but  not 
employing  any  corrections  like  those  of  the  last  formulation,  will  be 


presented.  To  this  end  we  introduce  the  quadratic  transverse  displacement 


wK  »  -  7  [x2  -  x2x  -  x3(x3  -  x2)  £  ,  y2  -  y^,  xy  -  x^]  <  (46) 

where 

s  *  §b  a  (47) 

and  Bb  Is  given  by  Eq.  (32).  In  the  corotational  description  adopted  here, 
the  displacement  w  vanishes  at  all  nodal  points.  The  approximation  of  the 
rotation  field  Is  linear,  as  It  Is  given  by  Eq.  (30). 

Upon  substltion  of  the  above  functions  into  Eqs.  (4)  and  (5),  one 
obtains  Bb  given  In  Eq.  (32)  and 

£  ■  5S  •  £  (W) 

where  B$  Is  given  by  Eq.  (33)  and 
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x3 
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Further  evaluation  of  the  above  equations  leads  to  the  conclusion  that 


which  proves  the  equivalence  between  the  previous  formulation  and  the  present 


one. 

The  structure  of  the  Eq.  (48)  clearly  shows  that,  within  this 
formulation,  ^  *  0  Indicates  that  the  Klrchhoff  mode  described  by  Eq.  (46)  is 
realized.  Thus,  in  the  thin-plate  limit,  quadratic  accuracy  could  be 
expected.  However,  all  the  remarks  made  with  regard  to  the  previous, 
equivalent  fqpnulation  are  also  pertinent  here.  This  means  that  for  some  mesh 
arrangements  the  results  obtained  with  the  present  formulation  are  very  poor 
unless  reduced  Integration  in  employed. 

V)  Mixed  formulation  with  the  modified  distribution  of  the  transverse 
displacements 

If  the  kinematics  presented  In  the  previous  formulation  and  the  constant 
shear  forces  related  to  the  matrix  Py  of  Eq.  (35)  are  used  in  Eq.  (8b) 
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where  if  Is  the  matrix  TT,  evaluated  at  the  centroid. 

The  above  result  Indicates  that  if  reduced  integration  is  to  be  applied 
along  with  the  modified  matrix  Bg,  the  centroid  should  be  selected  as  a  point 
of  Integration  rather  than  one  of  the  points  shown  on  Fig.  5(  in  [24]  the  word 
"centroid"  has  a  different  meaning). 
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The  Idea  here  is  the  same  as  in  the  related  section  concerning  beams. 
Namely,  for  given  nodal  rotations,  the  transverse  displacements  have  to  be 
found  which,  along  with  the  rotations,  form  the  bending  mode  that  is 
associated  with  no  shear  strain  energy.  This  displacement  mode  is  assumed  to 
be  defined  by  the  position  of  an  equivalent  Kirchhoff  configuration  in  which 
the  nodal  rotations  are,  in  an  average  sense,  closest  to  the  given 
rotations.  The  difference  between  the  given  nodal  displacements  and  the 
computed  ones  defines  the  shear  strain  energy  that  has  to  be  taken  into 
account. 

In  this  case  the  displacement  field  describing  basic  equivalent  Kirchhoff 
configuration  is  given  by  Eqs.  (46)  and  (47).  The  nodal  rotations  for  this 
field  are 


where  8^,  B^,  are  obtained  by  evaluating  the  matrix  B^  given  in  Eq.  (49) 
at  the  nodal  points  1,2,3  respectively.  Any  other  equivalent  Kirchhoff 
configuration  can  be  obtained  by  the  rigid  body  motion  given  in  Eq.  (46). 

Since  the  rigid  body  translation  is  here  Irrelevant  only  two  rotations  of  the 
magnitude  a„  and  a  will  be  considered.  In  this  case  the  difference  between 


given  nodal  rotations  and  those  of  an  equivalent  Klrchhoff  configuration  is 


A0  *  0*  +  a  dw  +  a  d„  -  0 
x  ~x  y  «y 


where 


dj  -  [1,  0,  1,  0,  1,  0] 

dj  -  [0,  1,  0,  1,  0,  1] 


(54a) 

(54b) 


The  function 


f(<V  <»y)  *  (40 )T  (40) 


Is  then  minimized  with  respect  to  ax,  ay.  The  result  is 


exi  -  8xi 


~y  3  '-yi  'yz  ^3'  '  ”yi  ey1  "  eyi 

The  quantities  on  the  right  hand  side  of  the  expressions  defining 
<*x  and  <iy  form  the  vector 

(0S)  *  C9xl.  0yr  0x2»  9y2»  ®x3’  0y33  ^ 


*1  *952  *9x3> 

■  9x1 

'll  ♦  “h  *  9?3> 

'  9y1 

(56a) 

(56b) 


which.  In  view  of  Eqs.  (30),  (33),  (48),  (50)  and  (52),  is 


The  matrices  IT.,,  1L,,  B,,  are  the  nodal  values  of  the  matrix  8„  of  Eq.  (45) 
'•SI  -St  -So  -s 

Because  of  the  corotational  description  adopted  here  the  initial  nodal 
displacements  are  zero,  so  the  transformation  of  the  bending  mode  to  the  total 
one  results  in 


y  m 


Cy) 


iay 


(59) 


So,  In  view  of  Eqs.  (55),  (57)  and  (59),  the  modified  matrix  §s  is 
T 


£s 


(60) 


and  It  describes  a  constant  distribution  of  the  shear  deformation  over  the 

element. 

Note,  that  matrix  f  s  of  Eq.  (45)  defines  a  linear  field  of  rotations 


over  the  element  and,  by  virtue  of  Eqs.  (56)  and  (58),  a  and  a  is  just  one 

x  y 


third  of  its  nodal  values.  Thus 


Therefore,  the  present  approach  once  again  confirms  that  in  [24]  the  centroid 
should  be  taken  as  the  integration  point.  More  importantly,  the  formulation 
seems  to  capture  the  predominant  mechanical  behavior  of  the  element  so  that  no 
reduced  integration  is  needed  in  conjunction  with  the  present  approach.  Its 
more  detailed  analysis  is  presented  in  [25]. 


5.  Numerical  examples 


Numerical  solutions  to  a  clamped  beam  and  to  a  simply  supported  circular 
plate,  both  under  uniform  loading,  are  presented.  If  C®  plate  elements  are 
used  in  the  analysis  of  thin  plates,  two  different  models  of  the  classical 
simply-supported  boundary  conditions  are  possible:  SSI,  In  which  only  the 
transverse  displacement  is  constrained  and  SS2,  in  which  both  the  transverse 
displacement  and  tangential  rotations  are  constrained.  In  the  present  case 
the  SSI  condition  has  been  employed.  Geometrical  and  mechanical  data  for  the 
beam  and  plate  problems  is  given  in  Tables  1  and  2,  respectively. 

The  central  displacements,  normalized  with  respect  to  the  analytic 
solution,  are  reported  in  Table  3  for  the  beam  and  Table  4  for  the  plate. 

The  number  of  elements  in  Table  3  refers  to  half  of  the  beam  since 
symmetry  is  used.  It  can  be  seen  that  the  Improvement  in  the  results  obtained 
by  changing  the  number  of  elements  from  1  to  2  Is  significant.  Doubling  the 
number  of  elements  to  4  does  not  change  the  displacement  much,  and  the 
accuracy  attained  Is  already  satisfactory.  The  difference  between  the  1- 
element  and  2-element  solutions  is  attributed  to  the  fact  that  a  single 
element  can  only  model  the  antisymmetric  mode  of  deformation  shown  in  Fig.  2, 
which  Is  associated  with  shear  strain  energy.  For  a  thin  structure,  this  is  a 
highly  energy-absorbing  mode  which  results  in  a  stiff  model. 

In  Table  4,  the  number  of  elements  is  for  a  quarter  of  the  plate;  the 
related  element  arrangements  are  shown  on  Fig.  6.  It  is  clear  that  the  first 
two  formulations  fall  while  formulation  (vi),  and  the  equivalent  formulation 
(v)  gives  the  best  results;  the  role  of  the  proper  decomposition  of  the  total 
C*^  configuration  into  Its  bending  and  shear  mode  is  therefore  apparent.  The 
formulations  (111)  and  (iv)  yield  results  equivalent  to  those  obtained  with 
one-point  quadrature  Introduced  in  [24].  However  a  different  selection  of  the 
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integration  point,  resulting  from  the  formulation  (v)  or  (vi 
significant  improvement,  especially  for  coarse  meshes. 


t ; 


I-. 


I 


),  gives 


-31- 


6.  Conclusions 

In  this  paper  two  basic  and  simple  C°  elements  have  been  investigated: 
the  beam  element  and  the  triangular  plate  element  (this  is  the  simplest  but  by 
no  means  the  easiest  plate  element);  both  with  linear  approximations  for 
rotations  and  displacements.  The  purpose  of  this  investigation  was  to  show 
that  the  proper  additive  decomposition  of  the  deformation  into  its  bending 
mode,  which  is  free  of  shear  strain  energy,  and  the  shear  mode,  is  crucial  for 
a  successful  development  of  C°  structural  elements.  This  is  clearly  seen  in 
the  analysis  of  the  triangular  linear  plate  element.  In  this  case,  almost  all 
the  formulations  are  different  and  in  most  cases  yield  unacceptable  results. 
Good  results  are  obtained  only  after  a  proper  definition  of  the  bending 
mode.  However  It  Is  Important  to  emphasize  that  the  proper  definition  of  the 
bending  mode  can  not  be  achieved  simply  through  the  use  of  reduced 
Integration. 

The  analysis  of  the  beam  element  shows  that  under  fortuitous 
circumstances  reduced  Integration  may  work.  In  the  beam  reduced  integration 
automatically  selects  the  proper  bending  mode  of  deformation.  This  Is 
probably  the  case  In  many  other  C°  elements,  as  for  instance  the  Lagrange 
family  of  plate  elements  which  are  based  on  reduced  integration.  Even  if  the 
optimal  bending  mode  of  deformation  Is  not  selected  by  reduced  integration,  in 
all  the  cases  where  this  approach  works,  the  reduced  integration  probably 
selects  a  bending  mode  which  Is  sufficiently  close  to  the  optimal  one  to  yield 
adequate  results.  There  are  however  cases, like  the  Serendipity  family  of 
plate  elements,  where  the  reduced  integration  fails,  [9,  12].  We  believe  that 
this  happens  because  of  inadequate  selection  of  the  bending  mode  of 
deformation. 

Of  the  formulations  presented  herein,  the  one  based  on  the  optimal 
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bending  configuration  and  originated  In  [25]  is  the  most  promising,  at  least 
for  simple  elements.  The  fact  that  it  leads  to  matrices  not  requiring  any 
reduced  integration  Indicates  that  it  removes  all  obstacles  to  the  correct 
behavior  of  the  elements  in  the  thln-structure  limit. 
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TABLE  1 


Data  for  the  beam  problem 

length 

height 

width 

Young  modulus 
Poisson  ratio 
shear  connection  factor 


10  in 
1  in 
1  in 

10.92  x  105  ps 
0.3 
5/6 


TABLE  2 

Data  for  the  plate  problem 


radius 

thickness 

Young  modulus 

Poisson  ratio 

shear  correction  factor 


5  in 
0.1  In 

10.92  x  105  ps 
0.3 
5/6 


TABLE  3 


h*. 


Central  displacement 
for  a  clamped  beam 


Notation: 

i)  reduced  integration  displacement  formulation 

ii)  mixed  formulation 

iii)  displacement  formulation  with  a  modified  matrix  B 

iv)  displacement  formulation  with  a  modified  distribution  of  the 
transverse  displacements 

v)  mixed  formulation  with  a  modified  distribution  of  the  transverse 
displacement 

vi)  displacement  formulation  based  on  the  optimal  bending  configuration 


TABLE  4 


Central  displacement 


for  simply  supported  circular  plate 


«*^No.  of  ele- 
— -~4ients 

Formul  atT5rr^-^__^ 

6 

24 

96 

0,  ii) 

0.063 

0.150 

0.398 

ill),  iv) 

0.722 

0.917 

0.981 

ref.  [24] 

0.703 

0.912 

0.948 

v),  vi) 

* 

0.824 

0.954 

0.989 

Notation  as  for  the  Table  3 


Fig.  1.  Sign  convention. 

Fig.  2.  Linear  beam  element: 

a)  total  deformed  configuration 
bl  symmetric  part  of  deformation 
c)  antisymmetric  part  of  deformation 

Fig.  3.  Linear  beam  element: 

a)  total  deformed  configuration 

b)  equivalent  Kirchhoff  configuration 

c)  shear  mode  of  deformation 

Fig.  4.  Geometry  of  triangular  plate  element. 

Fig.  5.  Position  of  the  Integration  point  in  ref.  [24] 
for  various  local  node  numbers. 

Fig.  6.  Discretization  of  the  circular  plate  example. 


Section  1 


INTRODUCTION 


In  Ref.  [1]  a  phenomena  called  membrane  locking  was  identified  and 
analyzed  for  curved  beam  elements  with  linear  axial  and  cubic  transverse 
displacement  fields.  Membrane  locking  results  from  the  inability  of  an 
element  to  bend  without  stretching:  since  any  bending  deformation  of  the 
element  is  then  accompanied  by  stretching  of  the  midline,  membrane  energy  is 
always  generated  in  bending,  thus  increasing  the  bending  stiffness.  Although 
the  study  employed  a  very  simple  beam  model,  it  is  clear  that  membrane  locking 
must  also  appear  in  shells,  (compare  [2-9]). 

However,  the  success  of  widely  employed  elements,  [10-15]  with  reduced 
shear  integration  and  the  claims  for  hybrid  and  mixed  elements  [16-19]  led  to 
the  conjecture  that  membrane  locking  may  be  circumvented  in  these  element:. 
Therefore,  a  similar  study  has  been  made  of  these  elements. 

It  has  been  found  that  reduced  shear  integration  in  curved  C°  elements 
does  mitigate  the  effects  of  membrane  locking.  Indeed,  a  complex 
interdependence  was  found  between  shear  and  membrane  underintegration: 
reduced  integration  in  either  the  shear  or  membrane  energies  leads  to  improved 
accuracy  in  the  bending  response.  In  fact  curved  elements  with  full 
membrane  integration  perform  quite  well  if  reduced  shear  integration  is 
used.  However,  reduced  shear  integration  is  accompanied  by  a  deterioration  of 
membrane-f lexural  coupling,  which  is  one  of  the  essential  features  of  a  curved 
element;  It  also  leads  to  the  appearance  of  kinematic  modes. 

Mixed  finite  elements  are  shown  to  also  exhibit  membrane  and  shear 
locking;  in  view  of  the  equivalence  ti.orems  [20], [21],  this  is  not 


surprising.  However,  for  a  beam,  hybrid  methods  are  not  subject  to  membrane 
locking  because  the  general  solution  to  the  equilibrium  equations  can  be 
constructed.  This  cannot  be  accomplished  for  arbitrary  shells,  but  it  may 
provide  the  insight  needed  for  a  rational  construction  cf  reduced-integration 
displacement  elements  which  avoid  locking.  We  are  convinced  such  displacement 
elements  provide  the  most  viable  approach  to  practical  computations;  the  extra 
calculations  associated  with  hybrid  and  mixed  elements  are  hard  to  justify 
when  displacement  elements  yield  the  same  results. 

In  Section  2,  the  governing  eauations  for  a  beam  based  on  shallow  shell 
theory  and  the  variational  formulations  which  pertain  in  this  context  to 
displacement,  hybrid,  and  mixed  elements  are  presented.  Using  a  specific  beam 
element,  the  interrelationship  of  membrane  and  shear  locking  is  demonstrated 
in  both  mixed  and  displacement  elements  in  Section  3.  Both  analytical  methods 
and  numerical  results  are  used.  In  Section  4,  isoparametric  beam  elements  are 
examined;  it  appears  that  the  cubic  element  avoids  locking,  though  the  cost  of 
its  counterpart  in  shell  analysis  is  quite  daunting.  Conclusions  are 
presented  in  Section  4. 


Section  2 


GOVERNING  EQUATIONS  AND  VARIATIONAL  FORMS 
Consider  a  curved  beam  of  height  h,  width  b  =  1,  which  is  approximated  by 
a  sequence  of  chords,  parametrized  by  x;  in  each  chord  the  shape  of  the  beam 
is  described  by  a  function  w(x)  as  shown  in  Fig.  1.  If  w(x)  is  small,  the 
behavior  of  the  beam  can  be  described  by  the  theory  of  shallow  structures 
wherein  the  following  equations  apply: 

i.  kinematic  relations  [22-25] 

e  =  u,x  +  w,x  v,x  (la) 

(lb) 

T  *  -♦  +  v.x  Uc) 

where  u  and  v  are  the  x  and  y  components  of  the  displacement  field  and  4  is 
the  rotation  of  the  cross-section;  z,  k,  and  y  are  the  membrane  (midplane) 
strain,  change  of  curvature  and  shear  deformation,  respectively;  commas  denote 
derivatives. 


constitutive  equations  (elastic) 

n  * 

Dj  *  Eh 

(2a) 

m  *  Dgic 

i)  .  ssi 
u2  12 

(2b) 

q  -  D3y 

*  icGh 

(2c) 

or 
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and  n,  m  and  q  are  the  membrane  (axial)  force,  moment  and  shear  on  the  cross- 
section  of  the  beam. 


iii.  equations  of  equilibrium 


"•x  ■  0 

• 

(3a) 

m,x  -  q  *  0 

(3b) 

»>x  +  <"w'x>"x  ’  P  ’ 

0 

(3c) 

boundary  conditions 

n  -  n\ 

or 

,  * 

(u  =  u  and 

5U  =  0) 

(4a) 

>x 

*E 

1 

II 

E 

or 

* 

( $  =  $  and 

6*  =  0) 

(4b) 

* 

q  +  nw,x  *  q  vx 

or 

/  * 

(v  =  v  and 

4v  =  0) 

(4c) 

where  asterisks  denote  prescribed  values,  the  prefix  &  a  variation,  and  vx  is 


4 


the  unit  normal  to  the  end,  which  takes  on  values  of  -1  and  +1  at  the  left  and 
right  hand  ends.  The  conditions  in  the  left  column  of  the  above  (which  are  on 
the  forces)  are  the  natural  boundary  conditions,  the  ones  on  the  right  the 
essential  boundary  conditions. 


Hu-Washizu  Functional 

The  Hu-Washizu  function  for  Eqs.  (1-3)  is  given  by 


2  *  *  + 

H  =  U-  fpvdx-£  (nu  +  m<*i  +  qv) 
C  i  =  l 


x  =  Xi 


[n  (e  -  u,x  -  w,xv,x)  +  m(K  +  $,x)  +  q(T 


(5a) 
v,x)]  dx 


where  Xj  and  ^  are  the  tw0  ends  0,:  the  beam  and  U  the  internal  energy,  which 
is  given  by 


U-I  /  (Die2  +  D2*2  +  D3y2)  dx 
2  C 


(5b) 


In  (5),  the  kinematic  variables  (u,  $,  v,  e,  k,  y)  are  all  independent,  as  are 
the  kinetic  variables  (n,  m,  q) ;  (u,  $,  v)  must  be  C°  functions; 

(e,  7*  n,  m,  q)  must  be  C”*,  C-*  functions  are  piecewise  continuous 

functions  which  are  allowed  to  be  discontinuous  across  element  interfaces. 

For  the  finite  element  formulation,  these  independent  variables  are 
approximated  by  shape  functions  S,  N,  and  E  as  follows 


5  £ 


(6) 


•  rt 
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(7) 


(8) 


where  d  is  the  matrix  of  nodal  displacements  and  a  and  e  are  the  discrete 
variables  for  the  stress  and  strain  fields,  respectively. 

We  also  introduce  the  strains  associated  with  the  displacement  field 
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where  the  kinematic  relations  (1)  are  used  in  the  first  equality  and  the 
elements  of  the  B  matrix  are  obtained  by  evaluating  the  expressions  in  the 
second  term  with  the  displacements  approximated  by  Eq.  (7). 

Substitution  of  Eqs.  (6-9)  into  (5)  provides  the  discrete  form  of  the  Hu 
Washizu  functional 
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This  functional  yields  the  following  discrete  equations: 


st rain -displacement 


B  d  *  E  e 


ii.  constitutive 


( lUe) 


(ion 


f T  r 


Iii.  equilibrium 


bt  6  -  fext 


Since  the  parameters  e  and  £  are  associated  with  C-1  functions,  the 

corresponding  element  matrices  ea  and  b  can  be  local  to  each  element. The 

^*6 

procedure  for  obtaining  the  governing  equations  is  then  the  following.  Using 
Eq.  (12),  It  follows  that 


*  r1  tT  U 


which,  when  introduced  into  Eq.  (11),  gives 
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B  =  (E  l)  E  )  B  d 


Now  the  use  of  equilibrium,  Eq.  (13)  leads  to  the  final  result 


K  d  =  fext 


k  •  y  lI  l 

L*  -w£ 

e 


(16a) 


(16b) 


where  L  is  the  convectivity  array,  the  element  stiffness  matrix  is  given  by 


k  -  bt  (i  r1  rT) _1  b 

■sg  <V  %<  -V 


The  heart  of  this  stiffness  lies  in  the  term  D,  see  Eq.  (10b),  which 
corresponds  in  form  to  the  usual  stiffness  except  that  E  has  replaced  B. 


Disolacement  Formulation  and  Selective  Reduced  Integration 


In  the  displacement  formulation 


k  I  (  m  }  *  D  /  « 


which  leads  to  the  fol lowing  identities  (compare  Eqs.  (6),  (8),  (9)) 


E  e  -  B  d 


S  B  =  D  B  d 

Si  S*  S* 


and  in  Eqs.  (10a)  the  coefficient  of  £  vanishes  and  e  can  be  eliminated 


■  *  - — r* 


give 


H  *  I  dT  K  d  -  dT  fext 

^  ■v  -w  •>*  •• 


(20a) 


and 


He 


D  B  dx 


(20b) 


This  stiffness  matrix  for  an  element  can  be  written  (see  Eqs.  (9)  and  (10g)) 


K  *  /  (D,  BT  B  ♦  09  B  +  D,  BT  B  )  dx 
~e  L  1  ~e  ~e  2  3  y 


membrane  flexural 


shear 


(21) 


Here  the  membrane,  flexural  and  shear  terms  are  identified.  The  crucial 
feature  of  a  curved  beam  is  that  w,x  does  not  vanish,  so  transverse 
displacements  v(x)  may  contribute  to  the  membrane  energy;  see  Eq.  (9). 
Moreover,  for  many  combinations  of  shape  functions,  any  transverse  reflection 
will  contribute  to  the  membrane  energy;  consequently  pure  bending 
deformations,  which  are  often  called  inextensional  modes  of  deformation,  can 
not  be  replicated  by  the  finite  element.  The  concomitant  increased  stiffness 
Is  called  "membrane"  locking  [1]. 

Mixed  Formulation 

In  the  mixed  formulation,  the  constitutive  equation  and  strain- 
displacement  equations  are  combined  so  that  the  independent  fields  are  the 
stresses  and  displacements.  Thus  Eqs.  (2d),  (6)  and  (8)  are  combined 
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(22) 


-1 


E  e  =  0  A  S  8 

"V  ^  ^ 


The  resulting  functional  H^j  is  obtained  by  substituting  Eq.  (22)  into  the 
first  and  last  term  of  Eq.  (10a),  which  gives 


\F  =  -  I  aT  IT1  6  +  fiT  F  d  -  dT  fext 

e  ^  t  ~  ~ 


Uc  =  f  ST  O'1  S  dx 


f  S'  0  1  S 
t  ~  ~  ' 


(23a) 


(23b) 


The  stationary  conditions  of  this  discrete  functional  are 


I  Se  '  2c  Se 

tf  Se  -  if* 

where 


F 


8  dx 


(24a) 

(24b) 


(24c) 


The  element  stiffness  is  obtained  by  combining  Eqs.  (24a)  and  (24b)  so  that  a 
relationship  is  obtained  between  ffi  and  d^  which  yields  that 


*e 


(25) 


The  possibility  of  locking  can  be  deduced  immediately  from  Eq.  (25).  Let 


10 


•Mr  MX  *i *5. - *  - 


'w 


S  =  D  B  (dim  8  »  dim  d)  (26) 

•W  -V  V 

Substituting  into  Eq.  (25),  we  obtain  a  K  which  is  identical  to  the 
displacement  formulation  stiffness,  Eq.  (20b).  This  result,  which  was 
obtained  in  a  more  sophisticated  way  in  [20],  shows  that  the  mixed  method  is 
equivalent  to  the  displacement  method  when  the  shape  functions  for  the 
stresses  are  obtained  by  Eq.  (26);  compare  the  conditions  for  equivalence  in 
[19].  Hence  locking  should  occur  in  mixed  models  (contrary  to  the 
implications  of  [19]),  and  reduced  integration  may  be  necessary;  this  will  be 
shown  later. 

The  counterpart  of  Eq.  (21)  ,  which  is  obtained  from  Eqs.  (24a)  and  (24b) 
is  that  for  a  mixed  model 

(27> 

Reduced  quadrature,  if  needed,  will  be  used  only  on  |;  all  terms  of  will 
be  integrated  exactly. 

Hybrid  Formulation 

In  the  hybrid  method,  all  stress  shape  functions  are  assumed  to  satisfy 
the  homogeneous  equilibrium  equations,  which  allows  the  second  term  in  Eq. 
(23a)  to  be  replaced  by  a  boundary  integral.  We  will  not  consider  this  form 
explicitly  but  instead  derive  the  hybrid  method  directly  from  Eq.  (23a). 
Details  are  given  in  the  next  Section. 


li 


Section  3 


CUBIC -LINEAR  ELEMENT 


Stiffness  Formulation. 

This  element,  shown  in  Fiy.  1,  employs  cubic  shape  functions  for  the 
transverse  deflection  v,  quadratic  shape  functions  for  the  rotations  of  the 
cross-section  $,  and  linear  shape  functions  for  the  axial  displacements  u. 
This  type  of  element  is  often  used  in  explicit  time  integration  because  its 
maximum  frequency  for  commonly  used  element  dimensions  is  lower  than  if  a 
cubic  is  used  for  u,  thus  providing  a  larger  stable  time  step  [9]. 

The  study  of  the  element  is  facilitated  by  considering  only  the  nodal 
degrees  of  freedom  which  are  associated  with  deformationai  modes,  thus 
excluding  all  rigid  body  modes.  The  deformationai  degrees  of  freedom  are 

d  *  CUgi  >  »  $2  *  »  <*2^  (28a) 

U21  =  u2  "  U1  (28b) 


where  U21  is  the  axial  displacement  of  the  right  end  relative  to  the  left 
one,  $.j  are  rotations  of  the  nodal  cross-sections  and  are  the  rotations  of 
the  tangents  to  the  middle  line  at  nodal  points.  The  shape  functions  and  the 
initial  shape  w(x)  are: 
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(2«d) 


w  =  L  a°  (c  -  2<;2  +  ;3)  +  L  (?3  -  £2) 

Note  that  these  shape  functions  exclude  rigid  body  motion  and  are  expressed  in 
a  corotational  system.  The  B  matrix  resulting  from  Eq.  (9),  is 
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(29) 

where  L  is  the  length  of  the  element  and  c  6  [u,l]  is  the  dimensionless 

parameter  of  its  chord.  The  orders  of  the  memDrane  and  shear  terms  for  this 
element  are  given  in  Table  1,  along  with  the  number  of  Gauss  quadrature  points 
required  for  exact  integration  of  these  polynomials. 

Analysis  of  Reduced  Integration  in  Displacement  Element 

Here  an  analysis  of  the  element  stiffness  will  be  presented  to 
demonstrate  the  effects  of  reduced  membrane  and  shear  integration.  Of  the 
rotational  degrees  of  freedom,  only  must  be  continuous  across  interelement 
boundaries,  so  is  local  to  each  element  and  can  be  eliminated  on  the 
element  level.  The  elimination  results  in  a  3x3  corotational  stiffness 

A 

matrix  K 


u 
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This  matrix  completely  defines  the  axial  and  flexural  properties  of  the 
element.  All  other  nodal  forces  are  found  through  global  equilibrium  of  the 


element.  The  entries  K ^ ,  K^3  >  »  *33  reflect  bending  properties  of  the 

element  while  represent  membrane -bending  coupling;  in  the  absence  of 

A  A 

this  coupling,  K12  =  K13  =  0.  To  make  the  analysis  tractable,  it  is  assumed 
that  Poission's  ratio  v  =  1/3  and 


0  _ 

-  a. 


=  -a(l+e) 


(31a) 


where  e  is  a  small  number.  Moreover,  we  assume  that 


0(1)  , 


a2  =  0(e) 


(31b) 


and  only  terms  of  order  e  are  retained  in  the  final  expressions.  These 
assumptions  are  satisfied  in  many  practical  applications 
where  -  0.1  and  a  -»  0.1.  Within  this  accuracy  we  obtain 

%1  “  *22  *1  +  *23  *2  ^32^ 

A  A 

The  coefficients  and  K23  >  for  various  reduced  integration  schemes  are 
given  in  Tables  2  and  3,  respectively. 

For  $2  *  0,  these  stiffness  terms  are  shown  in  Figs.  2  and  3  as  a 
function  of  the  element  slenderness  h/L.  As  expected,  for  a  fixed  value 
of  a,  the  bending  stiffness  of  the  curved  beam  approaches  that  of  the  straight 
beam  as  h/L  increases.  For  reduced  membrane  integration,  the  stiffness  is 


very  close  to  the  Euler-Bernoul 1 i  theory.  At  the  some  time,  a  clear 
difference  between  curved  and  straight  beams  is  retained. 

Reduced  shear  integration,  on  the  other  hand,  shifts  the  stiffness 
towards  that  of  a  straight  beam.  The  effects  of  curvature  appear  only  for 
very  slender  elements,  and  from  Figure  3,  it  is  apparent  that  this  is  also 
true  for  other  values  of  a.  Thus  reduced  integration  of  either  shear  or 
membrane  terms  reduces  the  bending  stiffness  and  yields  acceptable  results 
whenever  bending  is  the  predominant  mode. 

Remark  1.  For  the  case  o°  =  a^,  reduced  shear  integation  leads  to  singularity 
of  the  submatrix  which  has  to  be  inverted  to  eliminate  . 

Consequently,  cannot  be  uniquely  expressed  in  terms  of  the  remaining 
degrees  of  freedom.  The  nodal  forces  however,  in  this  case  depend  only  on  the 
mean  value  of  a.j  which  is  defined  uniquely.  To  avoid  this  difficulty  (which 
indicates  the  presence  of  a  zero  energy  mode)  this  submatrix  has  been 
integrated  exactly  in  this  analysis;  reduced  integration  is  used  for  all  other 
shear  contributions. 

Remark  2.  This  element  reduces  to  the  Euler-Bernoul 1 i  element  studied  in  Ref. 
[1]  if  Oj  =  <pr 


Hybrid  Formulation 

The  general  solution  of  the  homogeneous  form  of  Eqs.  (3)  is 


m  =  -0j  w  +  BgX  +  6^ 


The  parameters  g^,  g2,  S3  can  easily  be  interpreted  as  -n^,  -t^ ,  m^,  Fig. 
whereas  the  form  of  Eq.  (33)  results  in 


S  =  -w 


Note  that  the  development  of  a  general  equilibrium  solution  such  as  (33), 
while  feasible  for  curved  beams,  would  be  impossible  for  curved  shells. 

The  parameters  g^,  g2,  63  can  be  obtained  from  one  of  the  stationary 
conditions  of  the  functional  (23) ,  namely  Eq.  (24a)  ,  which  gives 


<£  sT  o-1  s  dx)  J  •  (£  sT  a)  dx  d 


The  right  hand  side  of  the  above  equation  can,  in  this  case,  be  evaluated 
without  assuming  any  approximation  for  the  displacement  field,  giving 


ST  Bdx)  d.^ST 


X-wk-w,xy 

A  A 

<x  +  y 


U2  -  U1 
v2  -  V1  -  *2L 
♦i  ‘  *2 


where  u^,  v^ ,  ^  are  translations  and  rotations  at  the  nodal  points.  Sine 
the  right  hand  side  of  Eq.  (35)  is  independent  of  the  approximation  of  the 


displacement  field,  we  may  claim  that  the  field  resulting  from  the  exact 
analytic  solution  of  the  problem  has  been  selected  and  that  the  element  will 
not  lock. 

The  equations  relating  the  nodal  forces  and  the  nodal  displacements  are 
obtained  by  solving  Eq.  (35)  for  ,  62*  63  and  then  making  use  of  the 
following  relationships 


(37) 


Mixed  Formulation 

For  the  cubic-linear  element,  the  following  internal  force  distribution 
were  considered: 


1.  nQ  -  qC;  this  distribution  conforms  with  the  guideline  of  Ref.  [19] 
that  the  internal  force  distribution  be  one  order  lower  than  that  obtained 
from  Eq.  (26): 

n  »  +  B2  U  -  5) 2  +  B3C2 

?  2 

m  «  b4L  (1  -  +  (38) 


*  H 


ii.  nQ  -  q L ;  in  tnis  element,  the  shear  distribution  is  consistent  in 
order  with  Eq.  ( 26) : 

n  -  Bj  +  S2  (1  -S)2  +  S352 

m  =  84  L  (1  -  5)  +  85  L  c  (39) 

q  ■  s6  +  67  (1  -  2?) 

iii.  nD  -  qC;  this  element  has  a  normal  force  distribution  consistent  with 
Eq.  (26): 

n  a  Bj  +  02  (-6c3  +  11c2  -  6c  +  1)  +  B3  (-6c3  +  7c2  -  2c) 
m  «  e4  L  (1  -  c)  +  85  L  C  (40) 

q  -  B6 

1v.  nC  -  qL;  this  element  has  the  simplest  normal  force  and  shear 
distribution  consistent  with  a  5  parameter  stress  model: 

n  - 

m  -  b2  L  (1  -  C)  +  S3  L  c  (41) 

q  *  S4  +  e5  (1  -  2c) 
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In  all  the  above  elements,  the  displacement  shape  functions  are  given  in 
Eq.  (28c)  and  the  B  matrix  needed  to  evaluate  Eq.  (24c)  is  given  in  Eq. 

(29).  In  reduced  membrane  integration,  the  rows  of  F  associated  witn  t‘  e 
axial  force  n  (e.g.  rows  I  to  3  for  nQ  -  qL)  are  undsrintegrated,  while  in  tne 

reduced  shear  integration,  the  rows  of  F  associated  with  tne  shear  q  (i.e. 

rows  6  and  7  for  nQ  -  qL)  are  underinteyrated. 

Numerical  Results 

Results  were  obtained  for  the  deep  arch  shown  in  Fig.  4.  The  results  for 

various  integration  schemes  in  displacement,  mixed  and  hybrid  formulations  are 

summarized  in  Table  4. 

The  displacement  method  with  full  integration  of  both  membrane  and  shear 
terms,  it  is  apparent,  is  far  too  stiff.  Reduced  integration  of  either  the 
membrane  or  shear  terms  leads  to  reasonable  agreement  with  the  analytic 
solution;  this  contrasts  with  the  behavior  of  curved  C1  elements  of  this  type, 
where  reduced  membrane  integration  Is  always  necessary  [1].  However,  reduced 
membrane  integration  with  full  shear  integration  gives  the  best  results.  The 
reason  for  this  is  apparent  from  examining  the  membrane-flexural  coupling 

A 

terms  K^,  which  vanish  for  shear  underintegration,  thus  providing  behavior 
similar  to  that  of  a  straight  element;  this  also  is  borne  out  by  comparing 
cases  3  and  5. 

The  hybrid  element  does  not  lock  and  is  in  good  agreement  with  the 
analytic  solution.  This  was  expected  in  view  of  the  fact  that  for  this 
element,  the  stiffness  matrix  is  Independent  of  the  shape  functions  used;  see 
comments  following  Eq.  (36). 

The  mixed  elements,  as  can  be  seen  from  Table  4,  also  exhibit  locking 
when  full  integrated  except  for  the  constant  membrane,  linear  shear  element 


(nC-qL).  In  the  quadratic  membrane,  linear  sheer  element  (nQ-qL) ,  various 
reduced  integration  schemes  were  tried.  The  best  accuracy  was  attained  for 
reduced  membrane  and  full  shear  integration.  Reduced  shear  integration  yields 
results  similar  to  that  of  the  quadratic  membrane,  constant  shear  (nQ-qC) 
curved  and  straight  elements,  supporting  the  hypothesis  that  reduce  shear 
integration  diminishes  flexural -membrane  coupling. 

This  conjecture  is  also  supported  by  the  effects  of  reduced  integration 

A 

on  the  stiffness  term  which  is  given  for  the  various  elements  in  Table 
4.  It  is  apparent  that  reduced  membrane  integration  provides  a  good  estimate 
of  this  stiffness  term,  while  for  reduced  shear  integration  it  often  vanishes, 
as  in  a  straight  beam. 

Examining  the  results  and  their  implications  in  more  detail  provides  some 
interesting  insight  into  the  relationship  between  reduced  integration  and 
shear/membrane  locking.  The  element  nQ-qC  (quadratic  n,  constant  q)  was 
designed  according  to  the  recommendation  of  [19];  see  Section  2.  As  can  be 
seen  from  Table  4,  nQ-qC  possesss  no  flexural -membrane  coupling  and  behaves 
like  a  straight  element  (compare  with  Case  13);  this  reflects  the  shear 
flexibility  brought  about  by  the  constant  shear  approximation.  The  element 
with  constant  n,  linear  q  (case  12)  in  fact  performs  better,  for  by  reducing 
the  order  of  n  and  increasing  the  order  of  q,  an  effect  similar  to  membrane 
underintegration  is  brought  about. 

These  results  (cases  7  and  12;  9  and  10)  again  illustrate  the 
Interrelationship  of  membrane  and  shear  underintegration  in  curved  elements, 
and  that  It  also  occurs  in  the  framework  of  mixed  elements. 

It  Is  worth  noting  that  when  continuous  distributions  were  used  for  the 
stress  functions  [18,  19],  no  locking  was  observed  despite  the  fact  that  the 
order  of  S  was  one  order  higher  than  that  given  by  Eq.  (26).  We  attribute 


CURVED  ISOPARAMETRIC  ELEMENTS 


Description  of  Elements 

The  element  studied  in  Section  3  employed  a  linear  axial  and  cubic 
transverse  displacement.  The  purpose  of  this  Section  is  to  examine  whether 
the  behavior  of  isoparametric  elements,  where  the  axial  and  transverse 
displacements  are  of  the  same  order,  is  similar. 

Two  elements  were  studied: 

1.  a  quadratic  isoparametric  with  3  nodes 

2.  a  cubic  isoparametric  with  4  nodes 

In  the  quadratic  and  cubic  elements,  u,  v,  $,  and  w  were  approximated  by 
quadratic  and  cubic  shape  functions,  respectively.  Because  the  curvature  is 
treated  by  shallow-shell  theory,  cf.  Eqs.  (1),  all  terms  were  integrated  over 
the  straight  x-axis.  Only  displacement  formulations  with  selective  reduced 
Integration  were  considered.  The  number  of  Gauss  points  required  for  exact 
quadrature  of  the  relevant  terms  is  given  in  Table  1. 

Some  of  the  reduced  quadrature  schemes  introduced  kinematic  modes; 
wherever  this  occured,  the  element  was  stabilized  by  exactly  integrating 
entries  of  the  stiffness  corresponding  to  the  transverse  displacements  of  the 
interior  nodes.  Although  this  scheme  is  not  suitable  for  practical 
applications,  it  provided  a  convenient  means  to  study  low  order  quadrature. 

All  results  obtained  in  this  manner  have  been  identified. 

Numerical  Results 

Results  for  the  curved  isoparametric  elements  are  presented  in  Tables  5 
and  6.  It  can  be  seen  from  Table  5  that  the  quadratic  element  exhibits  severe 


locking  when  full  quadrature  is  employed;  see  Table  1  for  the  quadrature  whicn 
is  exact.  Reduced  membrane  integration  alone  makes  a  big  difference  between  3 


and  2  Gauss  points  but,  still,  some  locking  is  present.  Further  reduction  to 
one  integration  point  practically  introduces  no  additional  change.  If  only 
shear  terms  are  under integrated ,  the  big  difference  is  between  2  and  1,  with 
almost  no  difference  between  3  and  2  Gauss  points. 

For  one-point  integration  of  the  shear  terms,  the  results  are  independent 
of  the  membrane  integration.  This  suggests  that  bending  of  the  element's 
midline  is  not  involved  in  the  deformation  process  and  the  membrane-bending 
coupling  is  almost  eliminated.  To  verify  this,  the  terms  K^.,  and  are 
given  in  Table  7  (see  Eq.  (30)).  It  is  clear  that  K^,  which  represents  the 
membrane -bending  coupling,  vanishes  for  one  point  shear  integration.  It  also 
vanishes  if  one  point  membrane  integration  is  used;  this  stems  from  the  fact 
that  the  initial  shape  of  the  element  is  symmetric  and  w,x  in  Eq.  (la) 
vanishes  at  the  center  of  the  element.  Obviously,  neither  one  point  shear  nor 
one  point  membrane  Integration  is  desirable. 

The  results  for  the  cubic  element  are  given  in  Table  6.  The  difference 
between  full  and  reduced  integration  is  quite  small.  This  element  is  only 
mildly  susceptible  to  membrane  or  shear  locking,  apparently,  deformation  modes 
which  eliminate  excessive  membrane  and  shear  energy  are  always  possible. 
However,  the  accuracy  is  best  for  3  point  quadrature  of  both  terms,  which 
represents  underintegration. 

Recommended  Quadrature  Scheme 

From  the  numerical  results,  it  appears  that  quadrature  schemes  which  in  a 
sense  filter  out  the  higher  order  terms  in  Eqs.  (la)  and  (lb),  are  most 
accurate.  For  an  isoparametric  element,  these  terms 


are  w,x  v,x  and  $,  respectively.  The  number  of  quadrature  points  would  then 
be  estimated  by  requiring  exact  quadrature  of  the  energies  associated  with  the 
lower  order  terms  in  Eqs.  (la)  and  (lb),  namely  u,x  and  v,x,  respectively. 

For  the  quadratic  isoparametric,  this  guideline  suggests  2  point 
quadrature  for  the  shear  and  membrane  terms,  while  for  the  cubic 
isoparametric,  it  suggest  3  points.  Tables  5  and  6  indicate  that  these 
quadrature  scheme  give  the  best  results. 

Remark  3.  Both  2x2  quadrature  in  the  quadratic  element  and  3x3  quadrature 
in  the  cubic  element  are  associated  with  kinematic  modes  in  the  three 
dimensional  shell  element.  These  modes  would  require  stabi 1 ization. 


CONCLUSIONS 


1.  Shear  membrane  locking  are  interrelated  in  curved  C°  elements.  Reduced 
integration  of  either  the  shear  or  membrane  terms  can  alleviate  locking, 
but  shear  underintegration  eliminates  the  membrane-f lexurai  coupling  which 
characterizes  curved  elements  and  thus  results  in  elements  whose 
performance  closely  approximates  that  of  straight  elements. 

2.  Mixed  finite  element  formulations,  when  used  with  generalized  stress 
fields  local  to  the  element,  also  exhibit  membrane  and  shear  locking. 

3.  In  cubic  isoparametric  beam  elements,  almost  no  locking  of  either  a 
membrane  or  shear  type  is  detectable  with  full  integration. 

4.  In  quadratic  isoparametric  beam  elements,  both  membrane  and  shear  locking 
are  present  and  the  interrelationship  between  these  types  of  locking 
described  in  the  first  conclusion  is  apparent. 

5.  Hybrid  curved  beam  elements  do  not  exhibit  locking.  This  is  probably  a 
consequence  of  the  fact  that  a  general  equilibrium  solution  can  be 
obtained  for  this  element;  it  is  doubtful  that  this  could  be  achieved  for 
a  curved  shell  element. 

It  Is  worth  noting  that  these  locking  phenomena  will  occur  regardless  of 


the  type  of  structural  or  continuum  theory  which  is  used.  A  shallow  beam 
theory  has  been  used  here  because  it  enables  the  order  of  the  membrane  terms 


to  be  easily  identified.  In  a  continuum  formulation,  these  terms  appear 
through  the  variations  in  the  Jacobian  which  are  Orouyht  about  by  the 
curvature,  and  are  not  as  easily  identified.  Nevertheless,  the  mechanical 
behavior  of  a  slightly  curved  element  will  be  identical.  Thus  in  either 
context,  the  use  of  higher  order  integration,  which  incidentally.,  is  often 
recommended  for  plastic  problems,  will  result  in  poor  element  performance 
because  of  locking. 

ACKNOWLEDGEMENT 

The  support  of  the  Air  Force  Office  of  Scientific  Research  under  Grant 
F4 9620 -82 -K -0013  under  the  direction  of  Dr.  A.  Amos  is  gratefully 
acknowledged. 


Fig.  1. 
Fig.  2. 


Notation  for  the  curved  beam  Element. 


Fig.  3. 
Fig.  4. 


Bending  stiffness  of  the  curved  beam  element  as  a  function  of 
aspect  ratio  normalized  with  respect  to  a  straight  beam;  F  and  R 
designate  full  and  reduced  quadrature  of  membrane  (M)  and  shear 
terms  (S);  EBA  designates  an  analytical  result  for  the  Euler- 
Bernoulli  beam. 

Bending  stiffness  of  the  curved  beam  element  as  a  functin  of  the 
initial  curvature  (see  Fig.  2  for  nomencl ature) . 

Problem  description  for  the  elastic  deep  aren. 
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TABLE  2 

Bending  stiffness  K ^  (Eq*  (30)) 

u  2  ,2 


(  4.0  +  Aj  a2  +  A2(^)  +  (^4  CA3  +  A4  a2  +  a5  (£)  ] 


Type  of 
integration 

A1 

a2 

*3 

A4 

i 

A5 

EBA 

0.0 

0.0 

0.5829 

0.0 

0.0 

FNI 

-0.8381 

-8.0 

0.9144 

0.8316 

1 .8288 

FS-RM 

-1.7984 

-2.6667 

0.6744 

0.0 

0.0 

FM-RS 

-0.8381 

' 

-8.0 

0.3143 

-0.0878 

-0.8381 

Notation:  EBA  -  Euler  Bernoulli  analytic  solution,  FNI  -  full  inte¬ 
gration,  FS-RM  -  full  shear  and  reduced  membrane  (2  Gauss 
points)  integration,  FM-RS  -  full  membrane  and  reduced 
shear  (1  Gauss  point)  integration. 


TABLE  3 


Bending  stiffness  <23  {Eq.  (30)) 


;  Eh 
'23  "  T2T 


|  2.0  +  Cj  a 


+  ^2  (•£ )  +  (-j^-)  [C3  +  +  Cs  (■ 


Type  of 
Integration 

ci 

C2 

C3 

C4  1 

t 

C5 
^  0 

EBA 

0.0 

0.0 

-0.2500 

0.0 

0.0 

FNI 

-0.8381 

-8.0 

-0.2856 

-0.8724 

l 

-0.8376 

FS-RM 

-1.7984 

-2.6667 

-0.6744 

0.0 

0.0 

FM-RS 

-0.8381 

-8.0 

0.3143 

-0.0878 

-0.8381 

Notation  the  same  as  for  Table  1. 
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TABLE  4 


Numerical  Results  for  a  Deep  Circular  Arch 
with  Cubic-Linear  Element 


<Case 

No 

Method 

Element  Used 

Further  Specifications 
of  the  Element 

(number  of  integration  points 
in  parenthesis) 

Force 

Displ 

lb 

in 

Stiffness 

terms** 

K22 

K12 

1 

Exact 

analytic 

EBA 

471.1 

— 

— 

2 

Displacement  Method 

8  curved 

shear 

flexible 

el. 

full  integration 

694.9 

1.25 

1.10 

3 

FM-RS,  full  mem.  (5)  red. 
shear  int. (1 ) 

487.8 

1.00 

0 

B 

FS-RM  full  shear,  (2) 
red.  mem.  int.  (2) 

473.3 

1.01 

1.12 

5 

8  straight 

EB  el. 

full  int. 

482.1 

0.9 

0.0 

6 

Hybrid  stress 

8  curved  shear 

473.3 

1.00 

1.00 

■ 

Mixed  Method 

8  el .  nQ  -  qC 

curved;  full  int. (4,1) 

484.5 

0.94 

0 

8 

8  el. 

nQ  -  qL 

curved;  full  int. (4,2) 

690.3 

1.19 

1.11 

9 

curved;  full  shear  (2) 
red.  mem.  (2) 

473.3 

1.01 

1.13 

10 

c 

curved;  full  mem.  (4) 
red.  shear  (1 )  • 

485.2 

0.94 

u 

•> 

8  el .  nD  -  qC 

curved;  full  int.  (4,1) 

487.0 

1.00 

0 

12 

8  el.  nC  -  qL 

curved;  full  int.  (3,2) 

473.4 

1.02 

1.13 

13 

8  el .  nC  -  qL 

straight;  ful 1  int. (1 ,2) 

482.1 

0.87 

0 

^Abbreviations:  EB  =Eu1er-Bernoulli ;  int.  =  integration;  mem.  =  membrane;  red.  *  reduced 
Normalized  with  respect  to  R22  and  <12  for  hybrid  stress  method  respectively 


TABLE  5 


Ratio  of  Numerical  to  Analytic  Results  for  a  Circular  Arch  Solution 
with  Quadratic  Isoparametric  Elements 
for  Reduced  Membrane  and  Shear  Quadrature 


(memb) 

Curved  Element 

Straight 

ng  (shearTx^^ 

3 

2 

1 

Element 

3 

1.575 

1.063 

1.084 

1.084 

2 

1.449 

1.004 

1  .023 

1.023 

1 

0.939* 

0.940* 

0.940* 

0.940* 

*)  Some  entries  of  the  stiffness  matrix  were  integrated  exactly  to 
stabilize  kinematic  modes. 


TABLE  6 

Ratio  of  Numerical  to  Analytic  Stiffness  fur  a  Deep  Circular  Arch  with 


Cubic  Isoparametric  Element  for  Reduced  Integration 


Gauss 

quadrature 
points  for 
membrane 

Gauss 

quadrature 
points  for 
shear 

Ratio 

5 

4 

1.010 

5 

3 

1.007 

5 

2 

1.004* 

5 

1 

0.957* 

4 

4 

1.010 

3 

4 

1 .004 

2 

4 

1.004 

1 

4 

1.023 

3 

3 

1.004 

*)  Kinematic  modes  were  stabilized. 

TABLE  7 

A  A 

Stiffnesses  and  K^  for  Quadratic  Isoparametric  Element 


'S'\^rig  (memb) 

nQ  (shear)\^ 

3 

2 

1 

CM 

CM 

*12 

*22 

*12 

*22 

*12 

3 

3.25 

1.09 

3.15 

1.12 

3.00 

0.00 

2 

1.11 

1.09 

1.01 

1.12 

0.87 

0.00 

1 

0.52 

0.00 

0.52 

0.00 

0.52 

0.00 

*)  normalized  with  respect  ot  hybrid  element  in  Table  3 
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